Light-meson leptonic decay rates in lattice QCD+QED

TL;DR

This paper presents the first lattice QCD+QED calculation of electromagnetic and isospin-breaking corrections to pion and kaon leptonic decay rates, refining CKM matrix element determinations and testing Standard Model unitarity.

Contribution

It introduces a novel lattice QCD+QED approach to compute decay rate corrections and improves the precision of CKM matrix element |V_{us}| and unitarity tests.

Findings

Electromagnetic corrections are 1.53(19)% for pions and 0.24(10)% for kaons.

The CKM matrix element |V_{us}| is determined as 0.22538(46).

First-row CKM unitarity is confirmed within uncertainties.

Abstract

The leading electromagnetic (e.m.) and strong isospin-breaking corrections to the $\pi^+ \to \mu^+ \nu[\gamma]$ and $K^+ \to \mu^+ \nu[\gamma]$ leptonic decay rates are evaluated for the first time on the lattice. The results are obtained using gauge ensembles produced by the European Twisted Mass Collaboration with $N_f = 2 + 1 + 1$ dynamical quarks. The relative leading-order e.m.~and strong isospin-breaking corrections to the decay rates are 1.53(19)\% for $\pi_{\mu 2}$ decays and 0.24(10)\% for $K_{\mu 2}$ decays. Using the experimental values of the $\pi_{\mu 2}$ and $K_{\mu 2}$ decay rates and updated lattice QCD results for the pion and kaon decay constants in isosymmetric QCD, we find that the Cabibbo-Kobayashi-Maskawa matrix element $ | V_{us}| = 0.22538(46)$, reducing by a factor of about $1.8$ the corresponding uncertainty in the Particle Data Group review. Our calculation of $|V_{us}|$ allows also an accurate determination of the first-row CKM unitarity relation $| V_{ud}|^2 + | V_{us}|^2 + | V_{ub}|^2 = 0.99988(46)$. Theoretical developments in this paper include a detailed discussion of how QCD can be defined in the full QCD+QED theory and an improved renormalisation procedure in which the bare lattice operators are renormalised non-perturbatively into the (modified) Regularization Independent Momentum subtraction scheme and subsequently matched perturbatively at $O(\alpha_{em}\alpha_s(M_W))$ into the W-regularisation scheme appropriate for these calculations.

Tables

Table 1. Table 1: Values of the coefficients Z η = 0 W-reg subscript superscript 𝑍 W-reg 𝜂 0 Z^{\textrm{W-reg}}_{\eta=0} (see Eq. ( 92 )) and Z fact superscript 𝑍 fact Z^{\mathrm{fact}} (see Eq. ( 93 )) calculated for the three values of the inverse coupling β 𝛽 \beta adopted in this work and for μ = 1 / a 𝜇 1 𝑎 \mu=1/a . In the fourth and sixth columns the values of the coefficient Z m fact superscript subscript 𝑍 𝑚 fact Z_{m}^{\mathrm{fact}} corresponding to the non-factorisable e.m. corrections to the mass RC in the MS ¯ ​ ( 2 ​ GeV ) ¯ MS 2 GeV \overline{\rm MS}(2\,\mbox{\rm GeV}) (see Eq. ( 40 )) are shown. The evaluation of the RCs in the RI'-MOM scheme has been carried out in Ref.  dicarlo using the methods M1 and M2 of Ref.  Carrasco:2014cwa (see Appendix  A ).

		Method M1		Method M2
$\beta$	$Z_{\eta =0}^{\text{W-reg}}$	$Z^{\mathrm{fact}}$	$Z_m^{\mathrm{fact}}$	$Z^{\mathrm{fact}}$	$Z_m^{\mathrm{fact}}$
$1.90$	0.00542 (11)	1.184 (11)	1.629 (41)	1.126 (7)	1.637 (14)
$1.95$	0.00519 (10)	1.172   (9)	1.514 (33)	1.123 (5)	1.585 (12)
$2.10$	0.00440   (7)	1.160   (6)	1.459 (17)	1.136 (4)	1.462   (6)

Table 2. Table 2: Values of the valence and sea bare quark masses (in lattice units), of the pion and kaon masses for the N f = 2 + 1 + 1 subscript 𝑁 𝑓 2 1 1 N_{f}=2+1+1 ETMC gauge ensembles used in Ref.  Carrasco:2014cwa and for the gauge ensemble, A40.40 added to improve the investigation of FVEs. A separation of 20 20 20 trajectories between each of the N cfg subscript 𝑁 cfg N_{\mathrm{cfg}} analysed configurations. The bare twisted masses μ σ subscript 𝜇 𝜎 \mu_{\sigma} and μ δ subscript 𝜇 𝛿 \mu_{\delta} describe the strange and charm sea doublet as in to Ref.  Frezzotti:2003xj . The values of the strange quark bare mass a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} , given for each β 𝛽 \beta , correspond to the physical strange quark mass m s phys ​ ( MS ¯ , 2 ​ GeV ) = 99.6 ​ ( 4.3 ) superscript subscript 𝑚 𝑠 phys ¯ MS 2 GeV 99.6 4.3 m_{s}^{\textrm{phys}}(\overline{\rm MS},2\,\mbox{{\rm GeV}})=99.6(4.3) MeV and to the mass RCs determined in Ref.  Carrasco:2014cwa . The central values and errors of pion and kaon masses are evaluated using the bootstrap procedure of Ref.  Carrasco:2014cwa .

ensemble	$\beta$	$V/a^4$	$N_{\mathrm{cfg}}$	$a{\mu }_{sea}=a{\mu }_{ud}$	$a{\mu }_{\sigma }$	$a{\mu }_{\delta }$	$a{\mu }_s$	$M_{\pi }\left(\mathrm{MeV}\right)$	$M_K\left(\mathrm{MeV}\right)$	$M_{\pi }L$
$A40.40$	$1.90$	${40}^3\times 80$	$100$	$0.0040$	$0.15$	$0.19$	$0.02363$	317 (12)	576 (22)	5.7
$A30.32$		${32}^3\times 64$	$150$	$0.0030$				275 (10)	568 (22)	3.9
$A40.32$			$100$	$0.0040$				316 (12)	578 (22)	4.5
$A50.32$			$150$	$0.0050$				350 (13)	586 (22)	5.0
$A40.24$		${24}^3\times 48$	$150$	$0.0040$				322 (13)	582 (23)	3.5
$A60.24$			$150$	$0.0060$				386 (15)	599 (23)	4.2
$A80.24$			$150$	$0.0080$				442 (17)	618 (14)	4.8
$A100.24$			$150$	$0.0100$				495 (19)	639 (24)	5.3
$A40.20$		${20}^3\times 48$	$150$	$0.0040$				330 (13)	586 (23)	3.0
$B25.32$	$1.95$	${32}^3\times 64$	$150$	$0.0025$	$0.135$	$0.170$	$0.02094$	259  (9)	546 (19)	3.4
$B35.32$			$150$	$0.0035$				302 (10)	555 (19)	4.0
$B55.32$			$150$	$0.0055$				375 (13)	578 (20)	5.0
$B75.32$			$80$	$0.0075$				436 (15)	599 (21)	5.8
$B85.24$		${24}^3\times 48$	$150$	$0.0085$				468 (16)	613 (21)	4.6
$D15.48$	$2.10$	${48}^3\times 96$	$100$	$0.0015$	$0.1200$	$0.1385$	$0.01612$	223  (6)	529 (14)	3.4
$D20.48$			$100$	$0.0020$				256  (7)	535 (14)	3.9
$D30.48$			$100$	$0.0030$				312  (8)	550 (14)	4.7

Equations

$$\Gamma (P^{\pm }\to {\ell }^{\pm }{\nu }_{\ell }[\gamma ])$$

$$\Gamma ({\pi }^{\pm }\to {\mu }^{\pm }{\nu }_{\ell }[\gamma ])$$

$$\Gamma (K^{\pm }\to {\mu }^{\pm }{\nu }_{\ell }[\gamma ])$$

$$S^{\text{full}}=\frac{1}{g_s^2}{S}^{\mathrm{YM}}+S^A+\sum _f\left\{S_f^{\text{kin}}+m_f{S}_f^m\right\}+\sum _{\ell }\left\{S_{\ell }^{\text{kin}}+m_{\ell }{S}_{\ell }^m\right\}.$$

$$R_1(aN;g_s,e,\text{m})=\frac{a{M}_{{\pi }^{+}}}{a{M}_{\Omega }}(aN;g_s,e,\text{m}),$$

$$R_2(aN;g_s,e,\text{m})=\frac{a{M}_{K^0}}{a{M}_{\Omega }}(aN;g_s,e,\text{m})$$

$$R_3(aN;g_s,e,\text{m})=\frac{a{M}_{D_s}}{a{M}_{\Omega }}(aN;g_s,e,\text{m}),$$

$$R_4(aN;g_s,e,\text{m})=\frac{a{M}_{K^{+}}-a{M}_{K^0}}{a{M}_{\Omega }}(aN;g_s,e,\text{m}),$$

$$R_5(aN;g_s,e,\text{m})=\frac{a{M}_{D^0}-a{M}_{D^{+}}}{a{M}_{\Omega }}(aN;g_s,e,\text{m}),$$

$$R_0(aN;g_s,e,\text{m})=\frac{a{M}_{\Omega }(aN;g_s,e,\text{m})}{M_{\Omega }^{\text{phys}}},$$

$$R_i(g_s,e,\text{m})\equiv \underset{N\to \infty }{\lim }{R}_i(aN;g_s,e,\text{m}),i=0,1,2,3,4,5.$$

$$R_i(g_s,e^{\text{phys}}(g_s),{\text{m}}^{\text{phys}}(g_s))=R_i^{\text{phys}},i=1,2,3,4,5.$$

$$a(g_s)=R_0(g_s,e^{\text{phys}}(g_s),{\text{m}}^{\text{phys}}(g_s)).$$

$$\frac{a{M}_H(L;g_s,e,\text{m})}{a{M}_H(g_s,e,\text{m})}=1-\kappa {\alpha }_{\mathrm{em}}{e}_H^2\left\{\frac{1}{2L{M}_H(g_s,e,\text{m})}+\frac{1}{L^2{M}_H^2(g_s,e,\text{m})}\right\},$$

$$S^{\text{QCD}}=\frac{1}{g_0^2}{S}^{\mathrm{YM}}+\sum _f\left\{S_{f,0}^{\text{kin}}+m_{f,0}{S}_f^m\right\},$$

$$a_0=\frac{⟨a_0{M}_{\Omega }{⟩}^{\text{QCD}}}{M_{\Omega }^{\text{phys}}}\text{and}a=\frac{⟨a{M}_{\Omega }{⟩}^{\text{full}}}{M_{\Omega }^{\text{phys}}}\equiv a_0(1+\delta a).$$

$$O^{\text{phys}}\equiv \frac{⟨aO{⟩}^{\text{full}}}{a}\text{and}{O}^{\text{QCD}}\equiv \frac{⟨a_{0}O{⟩}^{\text{QCD}}}{a_0}.$$

$$O^{\text{phys}}=O^{\text{QCD}}+\frac{⟨a_0\delta O{⟩}^{\text{QCD}}}{a_0}-\frac{\delta a}{a_0^2}⟨a_{0}O{⟩}^{\text{QCD}},$$

$$\frac{\delta M_{H_i}}{M_{H_i}}-\frac{\delta M_{\Omega }}{M_{\Omega }}=0(i=1,2,3,4).$$

$$\delta O^{\text{SIB}}=O^{\text{QCD}}-O^{\text{ISO}},$$

$$g(\mu )$$

$$m_f(\mu )$$

$$Z_m(g_0,a(g_s)\mu )m_{ud,0}(g_0)=\frac{{\hat{m}}_u(\mu )+{\hat{m}}_d(\mu )}{2}.$$

$$R_i^{\text{QCD-GRS}}=R_i^{\text{phys}}(1+{ϵ}_i^{\text{GRS}}),$$

$$\Gamma (P^{\pm }\to {\ell }^{\pm }{\bar{\nu }}_{\ell }[\gamma ])={\Gamma }^{\mathrm{QCD}}\cdot \left[1+{\overline{\delta R}}_P\right]+\mathcal{O}\left[{\alpha }_{\mathrm{em}}^2,(m_d-m_u{)}^2,{\alpha }_{\mathrm{em}}(m_d-m_u)\right],$$

$${\Gamma }^{\mathrm{QCD}}=\frac{G_F^2}{8\pi }|V_{q_1{q}_2}{|}^2{m}_{\ell }^2{\left(1-\frac{m_{\ell }^2}{M_P^{(0)2}}\right)}^2{f}_P^{(0)2}{M}_P^{(0)},$$

$$A_P^{(0)}\equiv ⟨0|{\bar{q}}_2{\gamma }_0{\gamma }_5{q}_1|P^{(0)}⟩\equiv f_P^{(0)}{M}_P^{(0)},$$

$${\cal H}_{W}=\frac{G_{F}}{\sqrt{2}}\,V_{q_{1}q_{2}}^{\ast}~{}O_{1}=\frac{G_{F}}{\sqrt{2}}\,V_{q_{1}q_{2}}^{\ast}\,\big{(}\bar{q_{2}}\gamma_{\mu}(1-\gamma_{5}){q_{1}}\big{)}\,\big{(}\bar{\nu}_{\ell}\gamma^{\mu}(1-\gamma_{5})\ell\big{)}\,,$$

$$O_{1}=\big{(}\bar{q_{2}}\gamma_{\mu}(1-\gamma_{5}){q_{1}}\big{)}\,\big{(}\bar{\nu}_{\ell}\gamma^{\mu}(1-\gamma_{5})\ell\big{)}\,,$$

$$\Gamma (P^{\pm }\to {\ell }^{\pm }{\bar{\nu }}_{\ell }[\gamma ])={\Gamma }^{(0)}\cdot \left[1+\delta R_P\right]+\mathcal{O}\left[{\alpha }_{\mathrm{em}}^2,(m_d-m_u{)}^2,{\alpha }_{\mathrm{em}}(m_d-m_u)\right],$$

$${\Gamma }^{(0)}=\frac{G_F^2}{8\pi }|V_{q_1{q}_2}{|}^2{m}_{\ell }^2{\left(1-\frac{m_{\ell }^2}{M_P^2}\right)}^2{\left[f_P^{(0)}\right]}^2{M}_P,$$

$${\mathcal{F}}_P^2\equiv \frac{\Gamma (P^{\pm }\to {\ell }^{\pm }{\bar{\nu }}_{\ell }[\gamma ])}{\frac{G_F^2}{8\pi }|V_{q_1{q}_2}{|}^2{m}_{\ell }^2{\left(1-\frac{m_{\ell }^2}{M_P^2}\right)}^2{M}_P}={\left[f_P^{(0)}\right]}^2(1+\delta R_P)$$

$$R_1(aN;g_s,e,\text{m})=\frac{a{M}_{{\pi }^0}}{a{\mathcal{F}}_{\pi }}(aN;g_s,e,\text{m}),$$

Full text

I Introduction

In flavour physics the determination of the elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix CKM , which contains just 4 parameters, from a wide range of weak processes represents a crucial test of the limits of the Standard Model (SM) of particle physics. Inconsistencies with theoretical expectations would indeed signal the existence of new physics beyond the SM and subsequently a detailed comparison of experimental measurements and theoretical predictions would provide a guide towards uncovering the underlying theory beyond the SM. For this to be possible non-perturbative hadronic effects need to be evaluated as precisely as possible and in this paper we report on progress in improving the precision of lattice computations of leptonic decay rates by including radiative corrections and strong isospin-breaking (IB) effects. A summary of our results has been presented in Ref. Giusti:2017dwk ; here we expand on the details of the calculation and include several improvements, most notably the renormalisation of the four-fermion weak operators in the combined QCD+QED theory (see Sec. IV). We also discuss in some detail how one might define the QCD component of the full (QCD+QED) theory (see Sec. II). Although such a separate definition of QCD is not required in order to obtain results computed in the full theory, it is necessary if one wishes to talk about radiative (and strong IB) “corrections” to results obtained in QCD. For this we need to specify what we mean by QCD.

The extraction of the CKM elements from experimental data requires an accurate knowledge of a number of hadronic quantities and the main goal of large-scale QCD simulations on the lattice is the ab initio evaluation of the nonperturbative QCD effects in physical processes. For several quantities relevant for flavour physics phenomenology, lattice QCD has recently reached the impressive level of precision of ${\cal{O}}(1\%)$ or even better . Important examples are the ratio $f_{K}/f_{\pi}$ of kaon and pion leptonic decay constants and the $K_{\ell 3}$ vector form factor $f_{+}(0)$ FLAG , which play the central role in the accurate determination of the CKM entries $|V_{us}/V_{ud}|$ and $|V_{us}|$, respectively. Such lattice computations are typically performed in the isospin symmetric limit of QCD, in which the up and down quarks are mass degenerate ($m_{u}=m_{d}$) and electromagnetic (e.m.) effects are neglected ($\alpha_{\mathrm{em}}=0$).

Isospin breaking effects arise because of radiative corrections and because $m_{u}\neq m_{d}$; the latter contributions are usually referred to as strong isospin breaking effects. Since both $\alpha_{\mathrm{em}}$ and $(m_{d}-m_{u})/\Lambda_{\mathrm{QCD}}$ are of ${\cal{O}}(1\%)$, IB effects need to be included in lattice simulations to make further progress in flavour physics phenomenology, beyond the currently impressive precision obtained in isosymmetric QCD.

Since the electric charges of the up and down quarks are different, the presence of electromagnetism itself induces a difference in their masses, in addition to any explicit difference in the bare masses input into the action being simulated. The separation of IB effects into strong and e.m. components therefore requires a convention. We discuss this in detail in Sec. II, where we propose and advocate the use of hadronic schemes, based on taking a set of hadronic quantities, such as particle masses, which are computed with excellent precision in lattice simulations, to define QCD in the presence of electromagnetism.

In recent years precise lattice results including e.m. and strong IB effects have been obtained for the hadron spectrum; in particular for the mass splittings between charged and neutral pseudoscalar (P) mesons and baryons (see for example, Refs. deDivitiis:2013xla ; Borsanyi:2014jba ). The QED effects were included in lattice QCD simulations using the following two methods:

• •

QED is added directly to the action and QED+QCD simulations are performed at a few values of the electric charge and the results extrapolated to the physical value of $\alpha_{\mathrm{em}}$ (see, e.g., Refs. Borsanyi:2014jba ; Boyle:2017gzv ; Hansen:2018zre );

• •

the lattice path-integral is expanded in powers of the two small parameters $\alpha_{\mathrm{em}}$ and $(m_{d}-m_{u})/\Lambda_{\mathrm{QCD}}$. This is the RM123 approach of Refs. deDivitiis:2011eh ; deDivitiis:2013xla ; Giusti:2017dmp which we follow in this paper.

In practice, for all the relevant phenomenological applications it is currently sufficient to work at first order in the small parameters $\alpha_{\mathrm{em}}$ and $(m_{d}-m_{u})/\Lambda_{\textrm{QCD}}$. The attractive feature of the RM123 method is that it allows one naturally to work at first order in isospin breaking, computing the coefficients of the two small parameters directly. Moreover, these coefficients can be determined from simulations of isosymmetric QCD.

The calculation of e.m. and strong IB effects in the hadron spectrum has a very significant simplification in that there are no infrared (IR) divergences. The same is not true when computing hadronic amplitudes, where e.m. IR divergences are present and only cancel in well-defined, measurable physical quantities by summing diagrams containing real and virtual photons BN37 . This is the case, for instance, for the leptonic $\pi_{\ell 2}$ and $K_{\ell 2}$ and semileptonic $K_{\ell 3}$ decay rates. The presence of IR divergences requires a new strategy beyond those developed for the calculation of IB effects in the hadron spectrum. Such a new strategy was proposed in Ref. Carrasco:2015xwa , where the determination of the inclusive decay rate of a charged P meson into either a final $\ell^{\pm}\nu_{\ell}$ pair or a final $\ell^{\pm}\nu_{\ell}\gamma$ state was addressed.

The e.m. corrections due to the exchange of a virtual photon and to the emission of a real one can be computed non-perturbatively, by numerical simulations, on a finite lattice with the corresponding uncertainties. The exchange of a virtual photon depends on the structure of the decaying meson, since all momentum modes are included, and the corresponding amplitude must therefore be computed non-perturbatively. On the other hand, the non-perturbative evaluation of the emission of a real photon is not strictly necessary Carrasco:2015xwa . Indeed, it is possible to compute the real emission amplitudes in perturbation theory by limiting the maximum energy of the emitted photon in the meson rest-frame, $\Delta E_{\gamma}$, to a value small enough so that the internal structure of the decaying meson is not resolved. The IR divergences in the non-perturbative calculation of the corrections due to the exchange of a virtual photon are cancelled by the corrections due to the real photon emission even when the latter is computed perturbatively, because of the universality of the IR behaviour of the theory (i.e., the IR divergences do not depend on the structure of the decaying hadron). Such a strategy, which requires an experimental cut on the energy of the real photon, makes the extraction of the relevant CKM element(s) cleaner.

In the intermediate steps of the calculation it is necessary to introduce an IR regulator. In order to work with quantities that are finite when the IR regulator is removed, the inclusive rate $\Gamma(P^{+}\to\ell^{+}\nu_{\ell}[\gamma])$ is written as Carrasco:2015xwa

[TABLE]

where the subscripts $0,1$ indicate the number of photons in the final state, while the superscript pt denotes the point-like approximation of the decaying meson and $\mu_{\gamma}$ is an IR regulator. In the first term on the r.h.s. of Eq. (1) the quantities $\Gamma_{0}(L)$ and $\Gamma_{0}^{\textrm{pt}}(L)$ are evaluated on the lattice. Both have the same IR divergences which therefore cancel in the difference. We use the lattice size $L$ as the intermediate IR regulator by working in the QEDL Hayakawa:2008an formulation of QED on a finite volume (for a recent review on QED simulations in a finite box see Ref. Patella:2017fgk ). The difference $\left[\Gamma_{0}-\Gamma_{0}^{\textrm{pt}}\right]$ is independent of the regulator as this is removed Lubicz:2016xro . As already pointed out, since all momentum modes contribute to it, $\Gamma_{0}(L)$ depends on the structure of the decaying meson and must be computed non-perturbatively. The numerical determination of $\Gamma_{0}(L)$ for several lattice spacings, physical volumes and quark masses is indeed the focus of the present study.

In the second term on the r.h.s. of Eq. (1) P is a point-like meson and both $\Gamma_{0}^{\textrm{pt}}(\mu_{\gamma})$ and $\Gamma_{1}^{\textrm{pt}}(\Delta E_{\gamma},\mu_{\gamma})$ can be calculated directly in infinite volume in perturbation theory, using a photon mass $\mu_{\gamma}$ as the IR regulator. Each term is IR divergent, but the sum is convergent BN37 and independent of the IR regulator. In Refs. Carrasco:2015xwa and Lubicz:2016xro the explicit perturbative calculations of $\left[\Gamma^{\textrm{pt}}_{0}+\Gamma^{\textrm{pt}}_{1}(\Delta E_{\gamma})\right]$ and $\Gamma_{0}^{\textrm{pt}}(L)$ have been performed with a small photon mass $\mu_{\gamma}$ or by using the finite volume respectively, as the IR cutoffs.

In Ref. Giusti:2017dwk we have calculated the e.m. and IB corrections to the ratio of $K_{\mu 2}$ and $\pi_{\mu 2}$ decay rates of charged pions and kaons into muons Giusti:2017dwk , using gauge ensembles generated by the European Twisted Mass Collaboration (ETMC) with $N_{f}=2+1+1$ dynamical quarks Baron:2010bv ; Baron:2011sf in the quenched QED (qQED) approximation in which the charges of the sea quarks are set to 0. The ratio is less sensitive to various sources of uncertainty than the IB corrections to $\pi_{\mu 2}$ and $K_{\mu 2}$ decay rates separately. In this paper, in addition to providing more details of the calculation than was possible in Ref. Giusti:2017dwk , we do evaluate the e.m. and strong IB corrections to the decay processes $\pi_{\mu 2}$ and $K_{\mu 2}$ separately. Since the corresponding experimental rates are fully inclusive in the real photon energy, structure-dependent (SD) contributions to the real photon emission should be included, however according to the Chiral Perturbation Theory (ChPT) predictions of Ref. Cirigliano:2007ga these SD contributions are negligible for both kaon and pion decays into muons. The same is not true to the same extent for decays into final-state electrons (see Ref. Carrasco:2015xwa ) and so in this paper we focus on decays into muons. The SD contributions to $\Gamma_{1}$ are being investigated in an ongoing dedicated lattice study of light and heavy P-meson leptonic decays.

An important improvement presented in this paper is in the renormalisation of the effective weak Hamiltonian. To exploit the matching of the effective theory to the Standard Model performed in Ref. Sirlin:1981ie it is particularly convenient to renormalise the weak Hamiltonian in the W-regularisation scheme. The renormalisation is performed in two steps. First of all, the lattice operators are renormalised non-perturbatively in the (modified) Regularization Independent Momentum subtraction (RI′-MOM) scheme at $O(\alpha_{\mathrm{em}})$ and to all orders in the strong coupling $\alpha_{s}$. Because of the breaking of chiral symmetry in the twisted mass formulation we have adopted in our study, this renormalisation includes the mixing with other four-fermion operators of different chirality. In the second step we perform the matching from the RI′-MOM scheme to the W-regularisation scheme perturbatively. By calculating and including the two-loop anomalous dimension at $O(\alpha_{\mathrm{em}}\alpha_{s})$, the residual truncation error of this matching is of $O(\alpha_{\mathrm{em}}\alpha_{s}(M_{W}))$, reduced from $O(\alpha_{\mathrm{em}}\alpha_{s}(1/a))$ in our earlier work Carrasco:2015xwa .

The main results of the calculation are presented in Sec. VI together with a detailed discussion of their implications. Here, we anticipate some key results: After extrapolation of the data to the physical pion mass, and to the continuum and infinite-volume limits, the isospin-breaking corrections to the leptonic decay rates can be written in the form:

[TABLE]

where $\Gamma^{(0)}$ is the leptonic decay rate at tree level in the Gasser-Rusetsky-Scimemi (GRS) scheme which is a particular definition of QCD Gasser:2003hk (see Sec. II.2.2 below). The corrections are about 1.5% for the pion decays and 0.2% for the kaon decay, in line with naïve expectations. Taking the experimental value of the rate for the $K_{\mu 2}$ decay, Eq. (3) together with $\Gamma^{(0)}(K^{\pm}\to\mu^{\pm}\nu_{\ell})$ obtained using the lattice determination of the kaon decay constant we obtain $|V_{us}|=0.22567(42)$, in agreement with the latest estimate $|V_{us}|=0.2253(7)$, recently updated by the Particle Data Group (PDG) PDG but with better precision. Alternatively, by taking the ratio of $K_{\mu 2}$ and $\pi_{\mu 2}$ decay rates and the updated value $|V_{ud}|=0.97420\,(21)$ from super-allowed nuclear beta decays Hardy:2016vhg , we obtain $|V_{us}|=0.22538(46)$. The unitarity of the first row of the CKM matrix is satisfied at the per-mille level; e.g. taking the value of $V_{us}$ from the ratio of decay rates and $|V_{ub}|=0.00413(49)$ PDG , we obtain $|V_{ud}|^{2}+|V_{us}|^{2}+|V_{ub}|^{2}=0.99988(46)$. See Sec. VI for a more detailed discussion of our results and their implications.

The plan for the remainder of this paper is as follows. A discussion of the relation between the “full” QCD+QED theory, including e.m. and strong IB effects, and isosymmetric QCD without electromagnetism is given in Sec. II. We discuss possible definitions of QCD in the full QCD+QED theory, and in particular we define and advocate hadronic schemes as well as the GRS scheme which is conventionally used Gasser:2003hk . In Sec. III we present the calculation of the relevant amplitudes using the RM123 approach. The renormalisation of the bare lattice operators necessary to obtain the effective weak Hamiltonian in the $W$-regularization scheme is performed in Sec. IV, while the subtraction of the universal IR-divergent finite volume effects (FVEs) is described in Sec. V. The lattice data for the e.m. and strong IB corrections to the leptonic decay rates of pions and kaons are extrapolated to the physical pion mass, to the continuum and infinite volume limits in Sec. VI. Finally, Sec. VII contains our conclusions. There are four appendices. The lattice framework and details of the simulation are presented in Appendix A. Appendix B contains a detailed discussion of the relation between observables in the full theory and in QCD, expanding on the material in Sec. II. An expanded discussion of the renormalisation of the effective weak Hamiltonian, including electromagnetic corrections, is presented in Appendices C, which contains a general discussion of the non-perturbative renormalisation in the RI'-MOM scheme and D in which issues specific to the twisted mass formulation are discussed.

II Defining QCD in the full theory (QCD+QED)

Before presenting the detailed description of our calculation of leptonic decay rates, we believe that it is useful to discuss the relation between the “full ” QCD+QED theory, that includes explicit e.m. and strong isospin breaking effects, and QCD without electromagnetism (denoted in the following as the full theory and QCD, respectively).

The action of the full theory can be schematically written as

[TABLE]

Here $g_{s}$ is the strong coupling constant, $S^{\textrm{YM}}$ is a discretisation of the gluon action, $S^{A}$ is the preferred discretisation of the Maxwell action of the photon, $S_{f}^{\textrm{kin}}$ is the kinetic term for the quark with flavour $f$, including the interaction with the gluon and photon fields, $m_{f}S_{f}^{m}=m_{f}\sum_{x}\overline{q}_{f}(x)q_{f}(x)$ is the mass term, $S^{\textrm{kin}}_{\ell}$ and $S^{m}_{\ell}$ are respectively the kinetic and mass terms for the lepton $\ell$ (for details see Appendix B). For fermion actions which break chiral symmetry, such as the Wilson action, a counterterm is needed to remove the critical mass and $m_{f}S_{f}^{m}$ has to be replaced with $m_{f}S_{f}^{m}+m_{f}^{\textrm{cr}}S_{f}^{\textrm{cr}}$. A mass counterterm is in principle needed also in the case of the lepton, but at leading order in $\alpha_{em}$ the lepton critical mass can be ignored.

At the level of precision to which we are currently working it is only the full theory, as defined in Eq. (4), which is expected to reproduce physical results and that is therefore unambiguous. Nevertheless, a frequently asked question is what is the difference between the results for a physical quantity computed in the full theory and in pure QCD, and how big are the strong isospin-breaking effects compared to the e.m. corrections. We particularly wish to underline that in order to properly formulate such questions it is necessary to carefully define what is meant by QCD. It is naturally to be expected that in QCD alone physical quantities will not be reproduced with a precision of better than ${\cal{O}}(\alpha_{\mathrm{em}})\simeq 1\%$ and this of course is the motivation for including QED. In order to define what is meant by QCD at this level of precision it is necessary to state the conditions which are used to determine the quark masses and the lattice spacing. The separation of the full theory into QCD and the rest is therefore prescription dependent.

In Ref. deDivitiis:2013xla the subtle issue of a precise definition of QCD has been discussed by using the scheme originally proposed in Ref. Gasser:2003hk , which we refer to as the GRS scheme and which has been widely used deDivitiis:2013xla ; Giusti:2017dmp ; Giusti:2017dwk . In the following and in Appendix B we present an extended and detailed discussion by introducing the hadronic schemes. Indeed, in light of the fact that hadron masses can nowadays be computed very precisely, we strongly suggest using hadronic schemes in future lattice calculations of QED radiative corrections. At the end of this section we discuss the connection with the GRS scheme that we had adopted at the time in which this calculation was started and that, for this reason, has been used in this work. A summary of the ideas discussed here has already been presented in Ref. Giusti:2018guw .

II.1 Renormalisation of the full theory

The main difference in the steps required to renormalise the full theory compared to the procedure in QCD is the presence of a massless photon and the corresponding finite-volume (FV) corrections which appear as inverse powers of $L$, where $L$ is the spatial extent of the lattice and the volume $V=L^{3}$. By contrast, in QCD for leptonic and semileptonic decays the FV corrections are exponentially small in the volume. In the discussion below, if necessary, we imagine that the chiral Ward identities have been imposed to determine the critical masses $m_{f}^{\textrm{cr}}$ Bochicchio:1985xa .

A possible strategy in principle is the following:

Fix the number of lattice points $N$, e.g. $T=2aN$ and $L=aN$, where $T$ and $L$ are the temporal and spatial extents of the lattice and the lattice spacing $a$ will be determined later. (The specific choice $T=2L$ is convenient for illustration but not necessary for the following argument.) 2. 2.

Using a four-flavour theory for illustration, we now need to determine the four physical bare quark masses, the bare electric charge and the lattice spacing. To this end we need to compute six quantities, e.g. the five dimensionless ratios111An alternative procedure to determine the bare electric charge would be the evaluation of the hadronic corrections to a leptonic observable.

[TABLE]

as well as a dimensionful quantity, e.g. the mass of the $\Omega$ baryon, computed in lattice units, from which the lattice spacing can be determined after extrapolation to the infinite volume limit (see below):

[TABLE]

where $M_{\Omega}^{\textrm{phys}}=1.672$ GeV is the physical value of the mass of the $\Omega$ baryon. For illustration we are considering the masses of QCD$+$QED stable pseudoscalar mesons in the numerators of the dimensionless ratios (5) and using $M_{\Omega}^{\mathrm{phys}}$ to determine the lattice spacing, but of course other quantities can be used instead. For example, in the four flavour theory that we are considering here one can in principle avoid potentially very noisy baryon observables by using one of the charmed mesons masses already considered above to set the scale. The choice of setting the scale with a charmed-meson observable could however, generate significant cutoff effects and reduce the sensitivity to the charm mass. In Eqs. (5) - (6) we have used $aN$ instead of $L$ to highlight that the infinite-volume limit should be taken at fixed lattice spacing (see Eq. (7) below). The quantity m represents the vector of bare quark masses $\textbf{m}\equiv\{m_{u},m_{d},m_{s},m_{c}\}$. Note that in the RM123 strategy, since one works at first order in $\alpha_{\mathrm{em}}$, it is not necessary to impose a renormalisation condition to fix the e.m. coupling deDivitiis:2013xla ; Giusti:2017dmp . In this case the electric charge can simply be fixed to the Thomson limit, i.e. $e=\sqrt{4\pi/137.036}$, and $R_{5}$ becomes a predictable quantity. For the remainder of this section, we assume that we are working to $O(\alpha_{\mathrm{em}})$ and only consider the four ratios $R_{i}$ (i=1,2,3,4) as well as $R_{0}$ when discussing the calibration of the lattices. Notice also that at first order in $\alpha_{\mathrm{em}}$ the $\pi^{0}$ cannot decay in two photons, so that it can also be used in the calibration procedure (see section III below). 3. 3.

Up to this point the procedure is the standard one used in QCD simulations. The difference here is in the FV effects which behave as inverse powers of $L$. We therefore envisage extrapolating the ratios $R_{i}$ to the infinite-volume limit:

[TABLE] 4. 4.

For a given discretisation and choice of $g_{s}$, the physical bare quark masses, $\textbf{m}^{\textrm{phys}}(g_{s})$, and the electric charge, $e^{\textrm{phys}}(g_{s})$, are defined by requiring that the five ratios $R_{1,2,3,4,5}$ take their physical values

[TABLE]

In practice, of course, depending on the specific choice of the ratios $R_{i}$, this will require some extrapolations of results obtained at different values of the bare quark masses and electric charge. 5. 5.

The lattice spacing $a$ at this value of $g_{s}$ can now be defined to be

[TABLE]

Note that with such a procedure the bare parameters and the lattice spacing $a$ do not depend on the lattice volume. 6. 6.

At first order in isospin breaking, i.e. ${\cal{O}}(\alpha_{\mathrm{em}},m_{d}-m_{u})$, the renormalisation of the lepton masses is performed perturbatively, by requiring that the on-shell masses correspond to the physical ones. If one wishes to go beyond first order, when hadronic effects first enter, then the physical lepton masses should be added to the quantities used in the non-perturbative calibration. The bare lepton masses, together with the other parameters, should be chosen such that, in addition to satisfying the conditions in Eq. (5), the lepton-lepton correlators decay in time as $e^{-m_{\ell}t}$, where $m_{\ell}$ is the physical mass of the lepton $\ell$.

In Eq. (7) we have taken the infinite-volume limit of the computed hadron masses. By working in the QED${}_{\textrm{L}}$ finite-volume formulation of QED, if for each hadron $H$ the FV corrections of order ${\cal{O}}(e^{2}/(M_{H}L)^{3},e^{4})$ can be neglected, then the extrapolation to the infinite-volume limit can be avoided by making use of the formula Hayakawa:2008an ; Borsanyi:2014jba (similar formulae also exist for other finite-volume formulations of the theory Lucini:2015hfa )

[TABLE]

where $e_{H}$ is the charge of the hadron $H$ and $\kappa=2.837297\,(1)$ is a known universal constant (independent of the structure of the hadron $H$). Equation (10) can be used to determine the infinite-volume mass of the hadron $H$ from the value measured on the finite-volume $L^{3}$, up to corrections of order of ${\cal{O}}(e^{2}/(m_{H}L)^{3},e^{4})$. (In any case, even if one wishes to study the behaviour with $L$ by performing simulations at different volumes, the subtraction of the universal $O(e^{2}/(M_{H}L))$ and $O(e^{2}/(M_{H}L)^{2})$ terms using Eq. (10) is a useful starting point; the residual leading behaviour of hadronic masses is then of $O(e^{2}/(M_{H}L)^{3},e^{4})$.)

II.2 Defining observables in QCD

The procedure discussed in section II.1 provides a full framework with which to perform lattice simulations of QCD together with isospin-breaking effects including radiative corrections. Nevertheless, one may wish to ask how different are the results for some physical quantities in the full theory (QCD+QED) and in QCD alone. We stress again that, under the assumption that isospin breaking effects are not negligible, QCD by itself is an unphysical theory and requires a definition. Different prescriptions are possible and, of course, lead to different results in QCD. In this section we propose and advocate hadronic schemes, based on the nonperturbative evaluation of a set of hadronic masses in lattice simulations and contrast this with schemes based on equating the renormalised strong coupling and masses in some renormalisation scheme and at a particular renormalisation scale which have been used previously.

We recall that the QCD action is given by

[TABLE]

where the kinetic term only includes the gluon links and the subscripts 0 indicate that the bare coupling and masses are different from those in the full theory of Eq. (4). Indeed the two theories have different dynamics that, in turn, generate a different pattern of ultraviolet divergences. The difference in the bare parameters of the two theories, for all schemes used to define QCD, can in fact be ascribed to the necessity of reabsorbing the different ultraviolet singularities. In what follows we present two different approaches to making the choice of the parameters $g_{0}$ and $m_{f,0}$. Explicit details of the lattice action, discretised using the Wilson formulation for the fermions for illustration, are presented in Appendix B.1.

II.2.1 Defining observables in QCD: hadronic schemes

In hadronic schemes we choose a value of $g_{0}$ and determine the bare quark masses $\textbf{m}_{0}^{\textrm{phys}}$ and the lattice spacing $a_{0}$ imposing the same conditions as for the full theory for the ratios $R_{0,\dots,4}$ evaluated at vanishing electric charge, i.e. following steps 1 - 5 in Sec. II.1 without imposing any constraint on the ratio $R_{5}$. We repeat that, for illustration we define the bare quark masses and lattice spacing using the five ratios $R_{i}$, but other hadronic quantities could be used instead, both in the full theory and in QCD. These parameters differ by terms of order ${\cal{O}}(\alpha_{\mathrm{em}})$ from those in the full theory. For this discussion, we make the natural and convenient choice $g_{0}=g_{s}$. (In order to make the perturbative expansion in Eq. (125) the difference $g_{s}-g_{0}$ should be less than ${\cal{O}}(\alpha_{\mathrm{em}})$.) With this choice, the lattice spacings in QCD ($a_{0}$) and in the full theory ($a$) are therefore given by

[TABLE]

To illustrate the procedure imagine that we wish to calculate an observable $O$ of mass dimension 1, for example the mass of a hadron which has not been used for the calibration. The generalisation to other cases is straightforward and presented in Appendix B. At a fixed value of $g_{s}=g_{0}$, we denote the best estimate of the observable $O$, which is the one obtained in the full theory, by $O^{\textrm{phys}}$, and that obtained in QCD as defined above by $O^{\textrm{QCD}}$:

[TABLE]

We define the difference of the two as being due to QED effects, $\delta O^{\textrm{QED}}\equiv O^{\textrm{phys}}-O^{\textrm{QCD}}$. There are 3 contributions to $\delta O^{\textrm{QED}}$:

The first contribution comes from the fact that the covariant derivatives in the kinetic terms in  (4) and Eq. (11) are different. This generates the diagrams in the correlation functions which contain the explicit exchange of virtual photons. 2. 2.

The second contribution comes from the fact that the bare quark masses appearing in Eq. (4) and Eq. (11) are different. The corresponding quark-mass counterterms must therefore be inserted into the correlation functions used to determine $O^{\textrm{phys}}$. We stress that the need to include quark-mass counterterms is generic and arises from the requirement that the conditions being used to determine the quark masses must be satisfied both in the full theory and in QCD (for the hadronic scheme being used for illustration we impose that the conditions in Eq. (8) are satisfied in both theories). 3. 3.

Finally we must account for the difference in the lattice spacings $\delta a=a-a_{0}$ in the full theory and QCD.

Combining these contributions we arrive at

[TABLE]

where we have combined the contributions to the correlation functions from the exchange of virtual photons and from the insertion of the mass counterterms into $\langle a_{0}\delta O\rangle^{\textrm{QCD}}$.

The detailed derivation of Eq. (14) is presented in Appendix B but some further comments may be helpful here. The first term on the right-hand side is one that can be calculated within QCD alone. It has a well defined continuum limit as does the sum of all the terms in Eq. (14). This term allows us to define what is the difference between QCD (defined as above) and the full theory in the hadronic scheme: $\delta O^{\mathrm{QED}}=O^{\textrm{phys}}-O^{\textrm{QCD}}$.

An important feature of the RM123 approach which we follow in the numerical study presented below, is that the ${\cal{O}}(\alpha_{\mathrm{em}})$ terms are computed explicitly and so we do not have to take the difference between numerical calculations performed in the full theory and in QCD. Each of the terms on the right-hand side of Eq. (14) is calculated directly. We now explain the procedure in some more detail by assuming that terms of order ${\cal{O}}(\alpha_{\mathrm{em}}^{2})$ are negligible (the extension to higher orders in $\alpha_{\mathrm{em}}$ is straightforward).

Correlation functions corresponding to diagrams with the exchange of a virtual photon and to the insertion of the mass counterterms are already of ${\cal{O}}(\alpha_{\mathrm{em}})$ and are calculated directly in QCD. The term proportional to the time separation in the correlation functions gives us the mass shift $\delta M_{H_{i}}$ ($i=1,2,3,4$) and $\delta M_{\Omega}$ for the five masses (or mass differences) in the ratios $R_{i}$ ($i=1,2,3,4$) in Eq. (5); 2. 2.

In the hadronic scheme being used for illustration, we impose the condition that the four ratios $R_{i}=m_{H_{i}}/m_{\Omega}$ are the same in QCD and in the full theory. This corresponds to requiring that

[TABLE]

The QED contribution to the left-hand side is different from zero (and also ultraviolet divergent) and we require the terms proportional to the counterterms to cancel this contribution. We therefore (in principle) scan the values of the four mass counterterms $\delta m_{f}=m_{f}-m_{f,0}$ ($f=u,d,s,c$) until the four conditions (15) are satisfied. Also in this case no subtraction of results obtained in the full theory and in QCD is necessary. 3. 3.

Finally we determine the difference $\delta a\equiv a-a_{0}$ in the lattice spacing. Having determined the bare masses using item 2, we can calculate the shift in the $\Omega$ mass, $\delta M_{\Omega}$ due to both QED and the mass counterterms and use Eq. (12). Since $a\delta M_{\Omega}$ is calculated directly, there is again no subtraction.

We have devoted a considerable discussion to the definition of the isospin-breaking effects due to electromagnetism, $\delta O^{\textrm{QED}}$. Having done this, the subsequent definition of the strong isospin breaking effects is straightforward. To do this however, we need to define the isosymmetric theory (labelled by “ISO”) by imposing appropriate conditions to determine the bare quark masses and the lattice spacing. Since $m_{u}=m_{d}$, in the $N_{f}=2+1+1$ theory we need to determine only three quark masses and hence we only need three conditions, e.g. we can use the ratios $R_{1,2,3}$ in Eq. (5) to determine the physical bare quark masses. For the determination of the lattice spacing we have two options. The simplest one is to work in a mass–independent scheme and set the lattice spacing in the isosymmetric theory, $a_{0}^{\textrm{ISO}}$, equal to the one of QCD with $m_{u}\neq m_{d}$, i.e. $a_{0}^{\textrm{ISO}}=a_{0}$. Notice that this choice is fully consistent with renormalisation because the ultraviolet divergences of the theories that we are considering do not depend on the quark masses. Note however, that they do depend instead on the electric charge. The other option is that we set the lattice spacing in the isosymmetric theory by using $R_{0}$ in Eq. (9). The difference between the two options is due to cutoff effects that disappear once the continuum limit is taken consistently. The strong isospin breaking correction $\delta O^{\textrm{SIB}}$ to the observable $O$ can now be defined by

[TABLE]

where $O^{\textrm{ISO}}=\frac{\langle a_{0}^{\textrm{ISO}}O\rangle^{\textrm{ISO}}}{a_{0}^{\textrm{ISO}}}$ is the value of the observable obtained in isosymmetric QCD. With these definitions we have the natural relation $O^{\textrm{phys}}=O^{\textrm{ISO}}+\delta O^{\textrm{QED}}+\delta O^{\textrm{SIB}}$. We underline however that $\delta O^{\textrm{SIB}}$ depends on the quantities used for calibration, both in 4-flavour QCD and in isosymmetric QCD.

II.2.2 Defining QCD: the GRS scheme

A different prescription, called the GRS scheme, was proposed in Ref. Gasser:2003hk to relate the bare quark masses and bare coupling of QCD ($m_{f,0}$ and $g_{0}$) to those in the full theory ($m_{f}$ and $g_{s}$). This prescription has been adopted in Refs. deDivitiis:2013xla ; Giusti:2017dmp ; Giusti:2017dwk . In the GRS approach, instead of determining the bare parameters of QCD by requiring that the chosen hadronic masses in QCD are equal to their physical values, one imposes that the renormalised parameters in a given short–distance scheme (e.g. the $\overline{\footnotesize\textrm{MS}}$ scheme) and at a given scale are equal in the full and QCD theories.

A consistent procedure is the following:

The full theory is renormalised by using a physical hadronic scheme as discussed in subsection  II.1. This means that for each chosen value of $g_{s}$ we know the corresponding physical value of the bare electric charge $e^{\textrm{phys}}(g_{s})$ and of the lattice spacing $a(g_{s})$. 2. 2.

The renormalisation constants (RCs) of the strong coupling constant and of the quark masses are computed in a short–distance mass–independent scheme both in the full theory and in the theory at vanishing electric charge. 3. 3.

In order to set the bare parameters of QCD at a given value of the lattice spacing we now chose a matching scale $\mu$ and impose that the renormalised strong coupling constant and the renormalised quark masses are the same as in the full theory. In practice we might want to simulate QCD at the same values of the lattice spacing used in the full theory simulations. In this case the matching conditions are

[TABLE]

where   $\widehat{}$   indicates quantities in the QCD+QED theory. Notice that quarks with the same electric charge have the same RC, e.g. $Z_{m_{u}}(e,g_{s},\mu)=Z_{m_{c}}(e,g_{s},\mu)$, and that the quark mass RC at vanishing electric charge is flavour independent, $Z_{m_{f}}(0,g_{0},\mu)=Z_{m}(g_{0},\mu)$. 4. 4.

In order to define isosymmetric QCD by using this approach, the bare up–down quark mass is determined from

[TABLE]

Some remarks are in order at this point. The GRS scheme is a short–distance matching procedure that can also be used to match the theories at unphysical values of the renormalised electric charge and/or quark masses with the physical theory.

By following the procedure outlined above one can perform lattice simulations of the full theory and of (isosymmetric) QCD at the same value of the lattice spacing but, consequently, at different values of the bare strong coupling constant. This is different from the strategy outlined in the previous subsection where, by using hadronic schemes, it was more natural to chose the same value of the bare strong coupling at the price of having two different lattice spacings. The absence of the lattice spacing counterterm (see Eq. (14) above) in the GRS scheme is compensated from the presence of the counterterm $(1/g_{0}^{2}-1/g_{s}^{2})S^{\mathrm{YM}}$ originating from the difference of the bare strong coupling constants in the two theories.

A remark of some practical relevance concerns the possibility of implementing hadronically the GRS scheme. To this end, note that in the GRS scheme the dimensionless ratios $R_{i}$ will not be equal to the corresponding physical values and the difference can be parametrized as follows

[TABLE]

where the $\epsilon_{i}^{\textrm{GRS}}$ are order ${\cal{O}}(\alpha_{\mathrm{em}})$ and depend on the chosen matching scheme and also on the chosen matching scale. Once the $\epsilon_{i}^{\textrm{GRS}}$ (and hence the $R_{i}^{\textrm{QCD-GRS}}$) are known, for example from a particularly accurate lattice simulation, then they can be used in other lattice computations. The bare quark masses are then determined by requiring that the $R_{i}$ in (isosymmetric) QCD reproduce $R_{i}^{\textrm{QCD-GRS}}$ as given by Eq.(19), and, at this stage, the GRS scheme can be considered to be a hadronic one as it is defined in terms of non-perturbatively computed quantities (in this case meson masses). We stress however that this requires prior knowledge of the $\epsilon_{i}^{\textrm{GRS}}$.

Of course other schemes are also possible. In general, the $\epsilon_{i}$ provide a unifying language to discuss the different schemes for the definition of (isosymmetric) QCD in the presence of electromagnetism; in physical hadronic schemes the $\epsilon_{i}=0$ while in the GRS and other schemes they are of order ${\cal{O}}(\alpha_{\mathrm{em}})$. For later use, we make the simple observation that two schemes can be considered to be equivalent in practice if the $\epsilon_{i}$ in the two schemes are equal within the precision of the computations.

Although the GRS scheme is perfectly legitimate, we advocate the use of physical hadronic schemes in future lattice calculations. For lattice simulations of physical quantities, a non-perturbative calibration of the lattice is necessary in general, but the renormalisation required for the GRS conditions in Eq. (17) is not generally necessary (except perhaps for the determination of the renormalised coupling and quark masses themselves). Now that hadronic masses are calculated with excellent precision in lattice simulations and their values are well known from experimental measurements, it is natural to use hadronic schemes. By contrast, the renormalised couplings and masses are derived quantities which are not measured directly in experiments. In spite of this, as explained above, at the time that our computation was started we chose to use the GRS scheme. Of course the physical results in the full theory do not depend on this choice.

III Evaluation of the amplitudes

At first order in $\alpha_{\mathrm{em}}$ and $(m_{d}-m_{u})/\Lambda_{\mathrm{QCD}}$ the inclusive decay rate (1) can be written as

[TABLE]

where $\Gamma^{\mathrm{QCD}}$ is the tree-level decay rate given by

[TABLE]

and $M_{P}^{(0)}$ and $f_{P}^{(0)}$ are the mass and decay constant of the charged P-meson mass defined in isosymmetric QCD in the chosen scheme.

The decay constant $f_{P}^{(0)}$ is defined in terms of the matrix element of the QCD axial current $A_{P}^{(0)}$ (in the continuum) as

[TABLE]

where the initial state meson $P^{(0)}$ is at rest. The decay rate is obtained from the insertion of the lowest-order effective Hamiltonian

[TABLE]

as depicted in the Feynman diagram of Fig. 1, where the decay of a charged kaon is shown as an example. At lowest order in $\alpha_{\mathrm{em}}$ the two full dots in the figure represent the two currents in the bare four-fermion operator

[TABLE]

whereas at order $\alpha_{\mathrm{em}}$ they will denote the insertion of the renormalised operator in the W-regularisation as defined in Sec. IV.

In order to compare our results for the e.m. and strong IB corrections to those obtained in Ref. Cirigliano:2011tm and adopted by the PDG PDG ; Rosner:2015wva however, we will use a modified expression:

[TABLE]

where $\Gamma^{(0)}$ is given by

[TABLE]

and $M_{P}$ is the physical mass of the charged P-meson including both e.m. and leading-order strong IB corrections.

The quantity $\delta R_{P}$ encodes both the e.m. and the strong IB leading-order corrections to the tree-level decay rate. Its value depends on the prescription used for the separation between the QED and QCD corrections, while the quantity

[TABLE]

is prescription independent Gasser:2010wz to all orders in both $\alpha_{\mathrm{em}}$ and $(m_{d}-m_{u})$.

The quantity $\mathcal{F}_{\pi}$ may be used to set the lattice scale instead of the $\Omega$ baryon mass. The physical value $\mathcal{F}_{\pi}^{\textrm{phys}}$ can be obtained by taking the experimental pion decay rate $\Gamma(\pi^{-}\to\mu^{-}\bar{\nu}_{\mu}[\gamma])=3.8408(7)\cdot 10^{7}$ s*-1* from the PDG PDG and the result for $|V_{ud}|=0.97420(21)$ determined accurately from super-allowed $\beta$-decays in Ref. Hardy:2016vhg . Consequently, one may replace $M_{\Omega}$ with $\mathcal{F}_{\pi}$ (as the denominator of the ratios $R_{1,\dots,4}$ in Eqs. (5)), $M_{\pi^{+}}$ with $M_{\pi^{0}}$ in the ratio $R_{1}$ (when working at leading order in $\alpha_{\mathrm{em}}$) and set the electron charge directly to its Thomson’s limit (instead of using the ratio $R_{5}$), namely

[TABLE]

Note that for the present study we were unable to use $M_{\Omega}$ to determine the lattice spacing because the corresponding baryon correlators were unavailable. The choice of using $\mathcal{F}_{\pi}$ instead to set the scale clearly prevents us from being able to predict the value of $|V_{ud}|$. This is one of the reasons why we advocate the use of hadronic schemes with hadron masses as experimental inputs for future lattice calculations. However, as already explained above, in this work we renormalise the QCD theory using the same set of hadronic inputs adopted in our quark-mass analysis in Ref. Carrasco:2014cwa , since we started the present calculations using the RM123 method on previously generated isosymmetric QCD gauge configurations from ETMC (see Appendix A). The bare parameters of these QCD gauge ensembles were fixed in Ref. Carrasco:2014cwa by using the hadronic scheme corresponding to $M_{\pi}^{(0),\textrm{FLAG}}=134.98$ MeV, $M_{K}^{(0),\textrm{FLAG}}=494.2(3)$ MeV and $f_{\pi}^{(0),\textrm{FLAG}}=130.41(20)$ MeV, while $M_{D_{s}}^{(0)}$ was chosen to be equal to the experimental $D_{s}^{+}$-meson mass, $M_{D_{s}^{+}}=1969.0(1.4)$ MeV PDG . Note that in the absence of QED radiative corrections $\mathcal{F}_{\pi}$ reduces to the conventional definition of the pion decay constant $f_{\pi}^{(0)}$. The superscript FLAG has been used because the chosen values of three out of the four hadronic inputs had been suggested in the previous editions of the FLAG review FLAG . For this reason we refer to the scheme defined from these inputs as the FLAG scheme.

We have calculated the same input parameters (28) used in the FLAG scheme also in the GRS scheme (corresponding to the $\overline{\footnotesize\textrm{MS}}$ scheme at $\mu=2$ GeV) obtaining222These values differ slightly from those obtained in Ref. Giusti:2017dmp , since we have now included the non-factorisable corrections of order ${\cal{O}}(\alpha_{\mathrm{em}}\alpha_{s}^{n})$ (with $n\geq 1$) to the mass renormalisation constant (see the coefficient $Z_{m}^{\mathrm{fact}}$ in Eq. (40) and in Table 1 below). We take the opportunity to update Eqs. (8), (10), (14) and (15) of Ref. Giusti:2017dmp with $\epsilon_{\pi^{0}}=0.01(4)$, $\epsilon_{K^{0}}=0.01(2)$, $\delta M_{D^{+}}+\delta M_{D^{0}}=1.7(1.0)$ MeV and $\delta M_{D_{S}^{+}}=2.3(4)$ MeV.: $M_{\pi}^{(0),\mathrm{GRS}}=135.0(2)$ MeV, $M_{K}^{(0),\mathrm{GRS}}=494.6(1)$ MeV, $M_{D_{s}}^{(0),\mathrm{GRS}}=1966.7(1.5)$ MeV and $f_{\pi}^{(0),\mathrm{GRS}}=130.65(12)$ MeV (see Eq. (111) in Sec. VI below). Therefore, the values of the inputs determined in the GRS scheme differ at most by $\sim 0.15\%$ from the corresponding values adopted in Ref. Carrasco:2014cwa for the isosymmetric QCD theory and the differences are at the level of our statistical precision. Thus, the result of our analysis of the scheme dependence can be summarized by the conclusion that the FLAG and GRS schemes can be considered to be equivalent at the current level of precision. Nevertheless, we have used the results of this analysis to estimate the systematic error on our final determinations of the isospin breaking corrections $\delta R_{P}$ induced by residual scheme uncertainties (see the discussion at the end of Sec. VI).

In light of this quantitative analysis, given the numerical equivalence of the two schemes at the current level of precision, in the rest of the paper we shall compare our results obtained in the GRS scheme with the results obtained by other groups using the FLAG scheme and we shall not use superscripts to distinguish between the two schemes.

The correction $\delta R_{P}$, defined in Eq. (25), is given by (see Ref. Carrasco:2015xwa )

[TABLE]

where

1. i)

the term containing $\mbox{log}(M^{2}_{Z}/M^{2}_{W})$ comes from the short-distance matching between the full theory (the Standard Model) and the effective theory in the $W$-regularisation Sirlin:1981ie ; 2. ii)

the quantity $\delta\Gamma_{P}^{\mathrm{(pt)}}(\Delta E_{\gamma})$ represents the ${\cal{O}}(\alpha_{\mathrm{em}})$ correction to the tree-level decay rate for a point-like meson (see Eq. (1)), which can be read off from Eq. (51) of Ref. Carrasco:2015xwa . The cut-off on the final-state photon’s energy, $\Delta E_{\gamma}$, must be sufficiently small for the point like-approximation to be valid; 3. iii)

$\delta A_{P}$ is the e.m. and strong IB correction to the decay amplitude $P\to\ell\nu$ with the corresponding correction to the amplitude with a point-like meson subtracted (this subtraction term is added back in the term $\delta\Gamma_{P}^{\mathrm{(pt)}}(\Delta E_{\gamma})$, see Eq. (1)). 4. iv)

$\delta M_{P}$ are the e.m. and strong IB corrections to the mass of the P-meson. The correction proportional to $2\,\delta M_{P}/M_{P}^{(0)}$ is present because of the definition of $f_{P}^{(0)}$ in terms of the amplitude and of the meson mass in Eq. (22).

Since we adopt the qQED approximation, which neglects the effects of the sea-quark electric charges, the calculation of $\delta A_{P}$ and $\delta M_{P}$ only requires the evaluation of the connected diagrams. These are shown in Figs. 1 - 5 for the case of $K_{\ell 2}$ decays. At ${\cal{O}}(\alpha_{\mathrm{em}})$ the diagram in Fig. 1 corresponds to the insertion of the operator renormalised in the W-renormalisation scheme.

In Eq. (29) $\delta A_{P}$ and $\delta M_{P}$ contain both the e.m. and the strong IB leading-order corrections

[TABLE]

where $\delta A_{P}^{W}$ is the e.m. correction from both the matching of the four-fermion lattice weak operator to the W-renormalisation scheme and from the mixing with several bare lattice four-fermion operators generated by the breaking of chiral symmetry with the twisted-mass fermion action which we are using. Both the matching and the mixing will be discussed and calculated in Sec. IV. As already pointed out, the renormalised operator, defined in the W-renormalisation scheme, is inserted in the diagram of Fig. 1. As for the diagrams of Figs. 2 - 5, which are already of order ${\cal{O}}(\alpha_{\mathrm{em}})$ and ${{\cal{O}}((m_{d}-m_{u})/\Lambda_{\textrm{QCD}})}$, it is sufficient to insert the weak current operator renormalised in QCD only.

In Eqs. (30) and (31) the quantity $\delta A_{P}^{\mathrm{SIB}}$ ($\delta M_{P}^{\mathrm{SIB}}$) represents the strong IB corrections proportional to $m_{d}-m_{u}$ and to the diagram of Fig. 4(b), 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. deDivitiis:2013xla ; deDivitiis:2011eh ). The term $\delta A_{P}^{J}$ ($\delta M_{P}^{J}$) is generated by the diagrams of Fig. 2(a-c), $\delta A_{P}^{T}$ ($\delta M_{P}^{T}$) by the diagrams of Fig. 2(d-e), $\delta A_{P}^{P}$ ($\delta M_{P}^{P}$) by the diagrams of Fig. 3(a-b) and $\delta A_{P}^{S}$ ($\delta M_{P}^{S}$) by the diagrams of Fig. 4(a-b). The term $\delta A_{P}^{\ell}$ corresponds to the exchange of a photon between the quarks and the final-state lepton and arises from the diagrams in Fig. 5(a-b). The term $\delta A_{P}^{\ell,\rm{self}}$ corresponds to the contribution to the amplitude from the lepton’s wave function renormalisation; it arises from the self-energy diagram of Fig. 5(c). The contribution of this term cancels out in the difference $\Gamma_{0}(L)-\Gamma_{0}^{\textrm{pt}}(L)$ and could be therefore omitted, as explained in the following section. The different insertions of the scalar density encode the strong IB effects together with the counter terms necessary to fix the masses of the quarks. The insertion of the pseudoscalar density is peculiar to twisted mass quarks and would be absent in standard Wilson (improved) formulations of QCD.

In the following subsection we discuss the calculation of all the diagrams that do not involve the photon attached to the charged lepton line. The determination of the contributions $\delta A_{P}^{\ell}$ and $\delta A_{P}^{\ell,\rm{self}}$ will be described later in subsection III.2.

III.1 Quark-quark photon exchange diagrams and scalar and pseudoscalar insertions

The terms $\delta A_{P}^{i}$ and $\delta M_{P}^{i}$ ($i=J,T,P,S$) can be extracted from the following correlators:

[TABLE]

where $\Delta_{\mu\nu}^{\mathrm{em}}(y_{1},y_{2})$ is the photon propagator, $J^{\rho}_{W}(x)$ is the local version of the hadronic $(V-A)$ weak current renormalised in QCD only333In our maximally twisted-mass setup, in which the Wilson $r$-parameters $r_{f_{1}}$ and $r_{f_{2}}$ are always chosen to be opposite $r_{f_{1}}=-r_{f_{2}}$ (see Appendix A), the vector (axial) weak current in the physical basis renormalises multiplicatively with the RC $Z_{A}$ ($Z_{V}$) of the axial (vector) current for Wilson-like fermions, i.e. $Z_{V}^{(0)}=Z_{A}$ and $Z_{A}^{(0)}=Z_{V}$ (see Appendix D).

[TABLE]

$j^{\mathrm{em}}_{\mu}$ is the (lattice) conserved e.m. current444The use of the conserved e.m. current guarantees the absence of additional contact terms in the product $j^{\mathrm{em}}_{\mu}(y_{1})j^{\mathrm{em}}_{\nu}(y_{2})$.

[TABLE]

and $T_{\mu}^{\mathrm{em}}$ is the tadpole operator

[TABLE]

In Eqs. (32) - (35) $\phi^{\dagger}_{P}(\vec{x},-t)=i\overline{q}_{f_{1}}(\vec{x},-t)\gamma_{5}q_{f_{2}}(\vec{x},-t)$ is the interpolating field for a P-meson composed by two valence quarks $f_{1}$ and $f_{2}$ with charges $e_{1}e$ and $e_{2}e$. The Wilson $r$-parameters $r_{f_{1}}$ and $r_{f_{2}}$ are always chosen to be opposite $r_{f_{1}}=-r_{f_{2}}$ (see Appendix A). We have also chosen to place the weak current at the origin and to create the P-meson at a negative time $-t$, where $t$ and $T-t$ are sufficiently large to suppress the contributions from heavier states and from the backward propagating P-meson (this latter condition may be convenient but is not necessary). In Eq. (35) $Z_{m}^{(0)}$ is the mass RC in pure QCD, which for our maximally twisted-mass setup is given by $Z_{m}^{(0)}=1/Z_{P}^{(0)}$, where $Z_{P}^{(0)}$ is the RC of the pseudoscalar density determined in Ref. Carrasco:2014cwa . The quantity ${\cal{Z}}^{f}_{m}$ is related to the e.m. correction to the mass RC

[TABLE]

and can be written in the form

[TABLE]

where ${\cal{Z}}_{\mathrm{QED}}^{f}$ is the pure QED contribution at leading order in $\alpha_{\mathrm{em}}$, given in the $\overline{\footnotesize\textrm{MS}}$ scheme at a renormalisation scale $\mu$ by Martinelli:1982mw ; Aoki:1998ar

[TABLE]

where $e_{f}$ is the fractional charge of the quark $q_{f}$ and $Z_{m}^{\mathrm{fact}}$ takes into account all the corrections of order ${\cal{O}}(\alpha_{s}^{n})$ with $n\geq 1$.

The quantity $Z_{m}^{\mathrm{fact}}$ is computed non-perturbatively in section IV and represents the QCD corrections to the “naive factorisation” approximation ${\cal{Z}}_{m}^{f}={\cal{Z}}_{\mathrm{QED}}^{f}$ (i.e. $Z_{m}^{\mathrm{fact}}=1$) introduced in Refs. Giusti:2017dmp ; Giusti:2017jof .

Analogously, the term $\left[\delta A_{P}\right]^{\mathrm{SIB}}$ and $\left[\delta M_{P}\right]^{\mathrm{SIB}}$ can be extracted from the correlator

[TABLE]

where, following the notation of Ref. Giusti:2017dmp , we indicate with $\widehat{m}_{f}$ and $m_{f}$ the renomalized masses of the quark with flavour $f$ in the full theory and in isosymmetric QCD only, respectively. We stress again that the separation between QCD and QED corrections is prescription dependent and in this work we adopt the GRS prescription of Refs. deDivitiis:2013xla ; Giusti:2017dmp ; Giusti:2017dwk , where

[TABLE]

Thus, in Eq. (42), the only relevant quark mass difference is $\widehat{m}_{d}-m_{ud}=-(\widehat{m}_{u}-m_{ud})$, whose value in the $(\overline{\rm MS},2\,\mbox{GeV})$ scheme was found to be equal to $1.19\,(9)$ MeV Giusti:2017dmp using as inputs the experimental values of the charged and neutral kaon masses.

Following Ref. deDivitiis:2013xla we form the ratio of $\delta C_{P}^{i}(t)$ with the corresponding tree-level correlator

[TABLE]

and at large time distances $t$ we obtain ($i=J,T,P,S,QCD$)

[TABLE]

where

[TABLE]

is the coupling of the interpolating field of the P-meson with its ground-state in isosymmetric QCD. The term proportional to $\delta M_{P}^{i}$ in the r.h.s. of Eq. (45) is related to the e.m. and strong IB corrections of the meson mass.

The function in the square brackets on the r.h.s. of Eq. (45) is an almost linear function of $t$. Thus, the correction to the P-meson mass, $\delta M_{P}^{i}$, can be extracted from the slope of the ratio $\delta C_{P}^{i}(t)/C_{P}^{(0)}(t)$ and the quantity $\delta[G_{P}^{i}A_{P}^{i}]$ from its intercept. As explained in Ref. Carrasco:2015xwa , in order to obtain the quantity $\delta A_{P}^{i}$ the correction $\delta G_{P}^{i}$ is separately determined by evaluating a correlator similar to those of Eqs. (32) - (35), in which the weak operator $J_{W}^{\rho}p_{P}^{\rho}/M_{P}$ is replaced by the P-meson interpolating field $\phi_{P}$.

For illustration, in Fig. 6 we show the ratios $C_{P}^{i}$ for the charged kaon ($P=K$) obtained from the ensemble D20.48 (see Appendix A). The top panel contains the ratio $\delta C_{K}^{\mathrm{SIB}}(t)/C_{K}^{(0)}(t)$, the ratio $\delta C_{K}^{J}(t)/C_{K}^{(0)}(t)$ is shown in the middle panel and the ratios $\delta C_{K}^{T}(t)/C_{K}^{(0)}(t)$ and $\delta C_{K}^{P}(t)/C_{K}^{(0)}(t)$ are presented in the bottom panel.

We find: i) the contributions $\delta C_{K}^{T}(t)/C_{P}^{(0)}(t)$ and $\delta C_{K}^{P}(t)/C_{P}^{(0)}(t)$ are separately large, but strongly correlated, since the tadpole insertion dominates the values of the e.m. shift of the critical mass $\delta m_{f}^{\mathrm{crit}}$ (see Ref. Giusti:2017dmp ). In the chiral limit they would cancel, but at finite masses the sum is small and linear in $t$. Because of the correlations it can nevertheless, be determined quite precisely (see the bottom right-hand plot of Fig. 6) where the sum is presented on an expanded scale. ii) the time dependence of the ratio $\delta C_{K}^{J}(t)/C_{K}^{(0)}(t)$ is almost linear in the time interval where the ground state is dominant.

III.2 Crossed diagrams and lepton self-energy

The evaluation of the diagrams 5(a) - (b), corresponding to the term $\delta A_{P}^{\ell}$ in Eq. (30), can be obtained by studying the correlator Carrasco:2015xwa

[TABLE]

where $S^{\ell}(0,x_{2})$ stands for the free twisted-mass propagator of the charged lepton. For the numerical analysis we have found it to be convenient to saturate the Dirac indices by inserting on the r.h.s. of Eq. (47) the factor $\left[\overline{v}(p_{\ell})\gamma_{\sigma}(1-\gamma_{5})u(p_{\nu})\right]$, which represents the lowest order “conjugate” leptonic ($V-A$) amplitude, and to sum over the lepton polarizations. In this way we are able to study the time behaviour of the single function $\delta C_{P}^{\ell}(t)$.

The corresponding correlator at lowest order ($O(\alpha_{\mathrm{em}}^{0})$) is

[TABLE]

In Eqs. (47) and (48) the contraction between the weak hadronic current $J_{W}^{\rho}(0)$ [see Eq. (36)] and its leptonic ($V-A$) counterpart gives rise to two terms corresponding to the product of either the temporal or spatial components of these two weak currents, which are odd and even under time reversal, respectively. Thus, on a lattice with finite time extension $T$, for $t\gg a$ and $(T-t)\gg a$ one has

[TABLE]

where $s_{0}=-1$, $s_{1,2,3}=1$ and

[TABLE]

is the relevant leptonic trace evaluated on the lattice using for the charged lepton the free twisted-mass propagator and for the neutrino the free Wilson propagator in the P-meson rest frame [${p_{P}^{\sigma}=(M_{P},\vec{0})}$].

Similarly, for the lowest-order correlator one has

[TABLE]

where $A_{P}^{(0)}$ is the renormalised axial amplitude evaluated on the lattice in isosymmetric QCD in the P-meson rest frame, namely

[TABLE]

The effect of the different signs of the backward-propagating signal in Eq. (49) can be removed by introducing the following new correlators:

[TABLE]

where

[TABLE]

Thus, the quantity $\delta A_{P}^{\ell}/A_{P}^{(0)}$ can be extracted from the plateau of the ratio $\delta\overline{C}_{P}^{\ell}(t)/\overline{C}_{P}^{\ell(0)}(t)$ at large time separations, viz.

[TABLE]

Note that the diagrams in Fig. 5(a) - (b) do not contribute to the electromagnetic corrections to the masses of the mesons and therefore the ratio (55) has no slope in $t$ in contrast to the ratios (45). Moreover, the explicit calculation of $X_{P}^{\ell,j}$ on the lattice is not required.

In terms of the lattice momenta $a\widetilde{p}_{\ell}$ and $a\overline{p}_{\ell}$, defined as

[TABLE]

the energy-momentum dispersion relations for the charged lepton and the neutrino in the P-meson rest frame are given by

[TABLE]

The 3-momentum of the final-state lepton $\vec{p}_{\ell}$ ($\vec{p}_{\nu}=-\vec{p}_{\ell}$) must be chosen to satisfy the equation

[TABLE]

Thus, for any given simulated P-meson mass $M_{P}^{(0)}$, the 3-momentum $\vec{p}_{\ell}=|\vec{p}_{\ell}|(1,1,1)$, is calculated from Eq. (60) and is injected on the lattice using non-periodic boundary conditions Bedaque:2004kc ; deDivitiis:2004kq for the lepton field. A simple calculation yields

[TABLE]

In Fig. 7 we show the correlators $C_{\pi}^{\mu(0)}(t)$, $\delta C_{\pi}^{\mu}(t)$, $\overline{C}_{\pi}^{\mu(0)}(t)$ and $\delta\overline{C}_{\pi}^{\mu}(t)$ for $\pi_{\mu 2}$ decays, multiplied by the ground-state exponential. These were obtained on the gauge ensembles A40.24 and D30.48 of Appendix A.

The subtraction of the backward signals, needed for extracting directly the quantity $\delta A_{P}^{\ell}$ given by Eq. (54), is beneficial also for extending the time region from which $\delta A_{P}^{\ell}$ (as well as the ratio $\delta A_{P}^{\ell}/A_{P}^{(0)}$) can be determined.

The quality of the signal for the ratio $\delta\overline{C}_{P}^{\mu}(t)/\overline{C}_{P}^{\mu(0)}(t)$ is illustrated in Fig. 8 for charged kaon and pion decays into muons for the case of the ensembles B55.32 and D30.48.

The calculation of the correction due to the diagram 5(c) is straightforward, since it is obtained by simply multiplying the lowest order amplitude, $A_{P}^{(0)}$, by the one-loop lepton self-energy evaluated on the lattice.

IV Renormalisation of the effective Hamiltonian and chirality mixing

In this section we provide the basic formalism to derive the e.m. corrections to the RCs non-perturbatively; further details of the calculation will be presented in a forthcoming publication dicarlo . This procedure relates the bare lattice operators to those in the RI'-MOM (and similar) renormalisation schemes up to order ${\cal{O}}(\alpha_{\mathrm{em}})$ and to all orders in $\alpha_{s}$. We also improve the precision of the matching of the weak operator $O_{1}$ (see Eq. (24)) renormalised in the RI'-MOM scheme to that in the W-regularisation by calculating the coefficient of the term proportional to $\alpha_{\mathrm{em}}\alpha_{s}\log(M_{W}^{2}/\mu^{2})$ in the matching coefficient. Using the two-loop anomalous dimension thus determined, we can evolve the operator to the renormalisation scale of $M_{W}$. Following this calculation the error due to renormalisation is reduced from order ${\cal{O}}(\alpha_{\mathrm{em}}\alpha_{s}(1/a))$ to order ${\cal{O}}(\alpha_{\mathrm{em}}\alpha_{s}(M_{W}))$.

The effective Hamiltonian, including the perturbative electroweak matching with the Standard Model Sirlin:1981ie , can be written in the form

[TABLE]

where the term proportional to the logarithm has been already included in Eq. (29) and $O_{1}^{\textrm{W-reg}}(M_{W})$ is the operator renormalised in the W-regularisation scheme, which is used to regularise the photon propagator. Since the W-boson mass is too large to be simulated on the lattice, a matching of the lattice weak operator $O_{1}$ to the W-regularisation scheme is necessary. In addition, for lattice formulations which break chiral symmetry, like the one used in the present study, the lattice weak operator $O_{1}$ mixes with other four-fermion operators of different chirality.

IV.1 The renormalised weak operator in the W-regularisation scheme

In order to obtain the operator renormalised in the W-regularisation scheme, we start by renormalising the lattice four-fermion operator $O_{1}$ defined in Eq. (24) in the RI'-MOM scheme Martinelli:1994ty , obtaining $O_{1}^{\text{{RI\textquotesingle}}}(\mu)$, and then perturbatively match the operator $O_{1}^{\text{{RI\textquotesingle}}}(\mu)$ to the one in the W-regularization Carrasco:2015xwa

[TABLE]

The coefficient $Z^{\textrm{W-{RI\textquotesingle}}}\left(M_{W}/\mu,\alpha_{s}(\mu),\alpha_{\mathrm{em}}\right)$ can be computed by first evolving the operator in the RI' scheme to the scale $M_{W}$ and then matching it to the corresponding operator in the W-scheme. The coefficient can therefore be written as the product of a matching coefficient and an evolution operator

[TABLE]

Below we will only consider terms of first order in $\alpha_{\mathrm{em}}$ and, therefore we will consistently neglect the running of $\alpha_{\mathrm{em}}$.

We note that the original bare lattice operators and $O^{\textrm{W-reg}}_{1}(M_{W})$ are gauge invariant, and thus the corresponding matching coefficients are gauge invariant. This is not the case for $O_{1}^{\text{{RI\textquotesingle}}}(\mu)$ that instead depends not only on the external states chosen to define the renormalisation conditions, but also on the gauge. Consequently the matching coefficient $Z^{\textrm{W-{RI\textquotesingle}}}\left(\frac{M_{W}}{\mu},\alpha_{s}(\mu),\alpha_{\mathrm{em}}\right)$ and the evolution operator $U^{\textrm{{RI\textquotesingle}}}\left(M_{W},{\mu},\alpha_{\mathrm{em}}\right)$ are in general gauge dependent. However, at the order of perturbation theory to which we are working, the evolution operator turns out to be both scheme and gauge independent.

In the following, we discuss in turn the matching coefficient, $Z^{\textrm{W-{RI\textquotesingle}}}\left(1,\alpha_{s}(M_{W}),\alpha_{\mathrm{em}}\right)$, the evolution operator $U^{\textrm{{RI\textquotesingle}}}\left(M_{W},{\mu},\alpha_{\mathrm{em}}\right)$, and the definition of the renormalised operator $O_{1}^{\mathrm{{RI\textquotesingle}}}(\mu)$, which will be obtained non-perturbatively.

a) The matching coefficient. At first order (one loop) in $\alpha_{\mathrm{em}}$

[TABLE]

where the strong interaction corrections for the RI'-MOM operator vanish, at this order, because of the Ward identities of the quark vector and axial vector currents appearing in the operator $O_{1}$ in the massless limit. We recall that we currently do not include terms of $O(\alpha_{s}(M_{W})\alpha_{\mathrm{em}})$ in the matching coefficient $Z^{\textrm{W-{RI\textquotesingle}}}$.

b) The evolution operator. The evolution operator $U^{\textrm{RI}^{\prime}}\left(M_{W},{\mu},\alpha_{\mathrm{em}}\right)$ is the solution of the renormalisation group equation

[TABLE]

where $U^{\textrm{{RI\textquotesingle}}}\left(M_{W},{\mu},\alpha_{\mathrm{em}}\right)$ satisfies the initial condition $U^{\textrm{{RI\textquotesingle}}}\left(M_{W},M_{W},\alpha_{\mathrm{em}}\right)=1$, $\gamma(\alpha_{s},\alpha_{\mathrm{em}})$ is, in general, the anomalous dimension matrix Buras:1993dy ; Ciuchini:1993vr , although in our particular case it is actually a number (and not a matrix), and $\beta(\alpha_{s},\alpha_{\mathrm{em}})$ is the QCD $\beta$-function:

[TABLE]

with

[TABLE]

where $N_{f}$ denotes the number of active flavours, and $N_{u}$ and $N_{d}$ denote the number of up-like and down-like active quarks respectively so that $N_{f}=N_{u}+N_{d}$. We may expand $\gamma(\alpha_{s},\alpha_{\mathrm{em}})$ in powers of the couplings as follows

[TABLE]

where $\gamma_{se}^{(1)}$ has been previously calculated in Ref. Brod:2008ss . In the case of the operator $O_{1}$ both $\gamma_{s}^{(0)}$ and $\gamma_{s}^{(1)}$ vanish whereas

[TABLE]

It can be demonstrated that, in addition to the leading anomalous dimension $\gamma_{e}^{(0)}$, $\gamma_{se}^{(1)}$ is also independent of the renormalisation scheme, thus in particular it is the same in RI' and in the W-regularisation schemes. It is then straightforward to derive $U^{\textrm{{RI\textquotesingle}}}\left(M_{W},{\mu},\alpha_{\mathrm{em}}\right)$

[TABLE]

Note that at this order the evolution operator is independent of the QCD $\beta$-function. This is a consequence of the fact that the QCD anomalous dimension vanishes for the operator $O_{1}$.

Combining Eqs. (63) - (65) and (71) we obtain the relation between the operator $O_{1}$ in the W-regularization scheme and the one in the RI'  scheme,

[TABLE]

which is valid at first order in $\alpha_{\mathrm{em}}$ and up to (and including) terms of ${\cal{O}}(\alpha_{\mathrm{em}}\alpha_{s}(M_{W}))$ in the strong coupling constant.

c) The renormalised operator in the RI'-MOM scheme. When we include QCD and e.m. corrections at ${\cal{O}}(\alpha_{\mathrm{em}})$, the operator $O_{1}$ on the lattice with Wilson fermions mixes with a complete basis of operators with different chiralities. In addition to $O_{1}$, the mixing involves the following operators

[TABLE]

The mixing is a consequence of the explicit chiral symmetry breaking of Wilson-like fermions on the lattice. Therefore, the renormalised operators in the RI'-MOM scheme, $\vec{O}^{\textrm{{RI\textquotesingle}}}(\mu)$, with $\vec{O}=(O_{1},\ldots O_{5})$, can be written in terms of bare lattice operators $\vec{O}^{\mathrm{bare}}(a)$ as

[TABLE]

where $Z_{O}(\mu a)$ is a $5\times 5$ renormalisation matrix. We note that in pure QCD the operator $O_{1}$ mixes only with $O_{2}$, with scale independent coefficients, whereas the full $5\times 5$ renormalisation matrix is necessary in general when e.m. corrections are included.

We find it particularly convenient to rewrite Eq. (74) in the form

[TABLE]

where $Z^{\mathrm{QCD}}$ is the mixing matrix in pure QCD (corresponding to $\alpha_{\mathrm{em}}=0$), and

[TABLE]

is the pure, perturbative QED mixing matrix (corresponding to $\alpha_{s}=0$). In Eq. (75) we have introduced the ratio

[TABLE]

so that, at first order in $\alpha_{\mathrm{em}}$, Eq. (75) is written as

[TABLE]

The ratio ${\cal R}$ encodes all the nonperturbative contributions of order $\mathcal{O}(\alpha_{\mathrm{em}}\alpha_{s}^{n})$ with $n\geq 1$, other than the factorisable terms given by the product $Z^{\mathrm{QED}}\,Z^{\mathrm{QCD}}$. In other words if $Z_{O}$ were simply given by $Z_{O}=Z^{\mathrm{QED}}\,Z^{\mathrm{QCD}}$ at first order in $\alpha_{\mathrm{em}}$ then $\eta$ would be zero. The case $\eta=0$ thus corresponds to the factorisation approximation that was first introduced in Refs. Giusti:2017dmp ; Giusti:2017jof .

In this work, the ratio ${\cal R}$ has been computed non-perturbatively on the lattice to all orders in $\alpha_{s}$ and up to first order in $\alpha_{\mathrm{em}}$. Introducing this ratio ${\cal{R}}$ in the non-perturbative calculation is useful since by using the same photon fields in the lattice calculation of $Z_{O}$ and $Z^{\mathrm{QED}}$ the statistical uncertainty due to the sampling of the photon field is significantly reduced. Note that the ratio is also free from cut-off effects of ${\cal{O}}(\alpha_{\mathrm{em}}a^{n})$. The non-perturbative calculation of ${\cal{R}}$, in terms of the matrix $\eta$, is described in Appendix C and all the details and results will be presented in a forthcoming publication dicarlo .

As already mentioned, pure QCD corrections in Eq. (78) only induce the mixing of the operator $O_{1}$ with the operator $O_{2}$. This mixing produces the renormalised QCD operators

[TABLE]

which, similarly to the corresponding continuum operators, belong respectively to the $(8,1)$ and $(1,8)$ chiral representations with respect to a rotation of the quark fields Bochicchio:1985xa . These are the combinations entering on the r.h.s. of Eq. (78).

When we include the e.m. corrections at ${\cal{O}}(\alpha_{\mathrm{em}})$, the matrices $\Delta Z^{\mathrm{QED}}$ and $\eta$ in Eq. (78) induce, in general, the mixing of $O_{1}^{\chi}$ with the full basis of operators in Eq. (73). As shown in Appendix D, however, in the twisted-mass formulation used in this paper the only relevant chirality mixing is the one between the operators $O_{1}$ with $O_{2}$. Indeed, the mixing coefficients with the operators $O_{3}$ and $O_{4}$ are found to be odd in the parameter $\overline{r}\equiv r_{1}r_{\ell}=-r_{2}r_{\ell}$, defined by the product of the Wilson $r$-parameters of the valence quarks and the lepton (with $r_{2}=-r_{1}$ in our procedure). Therefore, taking the average over the values of the parameter $\overline{r}$ (with $\overline{r}=\pm 1$) when computing the amplitude, eliminates the mixing with $O_{3}$ and $O_{4}$. Moreover, the matrix element of the operator $O_{5}$ between a pseudoscalar meson and the vacuum vanishes, so that the mixing with the operator $O_{5}$ cannot contribute to the decay rate. Therefore, Eq. (78) for the renormalised operator $O_{1}^{\textrm{{RI\textquotesingle}}}$ simplifies to

[TABLE]

where we have explicitly indicated the dependence of the various terms on $\alpha_{s}$ and the renormalisation scale. Since the mixing of the bona fide $(8,1)$ operator $O_{1}^{\chi}$ with $O_{2}^{\chi}$ is a consequence of the explicit chiral symmetry breaking of Wilson-like fermions on the lattice, the corresponding coefficient is due to lattice artefacts and can only be a function of the lattice bare coupling constant $\alpha_{s}(1/a)$ Bochicchio:1985xa .

IV.2 Complete expression for the matching coefficients

We are now in a position to collect the results of the previous subsection in order to provide the final expression relating the renormalised operator $O^{\textrm{W-reg}}_{1}$ in the W-regularization to the lattice bare operators $O_{1}$ and $O_{2}$ at first order in $\alpha_{\mathrm{em}}$. Combining Eqs. (72) and (80) and choosing $\mu=1/a$ as renormalisation scale in the intermediate RI'-MOM scheme we obtain:

[TABLE]

Using the results of Ref. Carrasco:2015xwa , obtained in perturbation theory at order ${\cal{O}}(\alpha_{s}^{0})$, we have determined the values for the matching and mixing coefficients,

[TABLE]

where $\xi$ is the photon gauge parameter ($\xi=0\,(1)$ in the Feynman (Landau) gauge). It is worth noting that the renormalised operator in the W-regularization scheme is gauge independent, at any order of perturbation theory. In particular, as shown by Eq. (82), at first order in $\alpha_{\mathrm{em}}$ and at zero order in $\alpha_{s}$ the gauge dependence of the matching coefficient of $O_{1}^{\chi}$ cancels in the sum $C^{\textrm{W-{RI\textquotesingle}}}+\Delta Z^{\mathrm{QED}}_{11}=-15.539$. By contrast, for the matching coefficient of $O_{2}^{\chi}$, the two terms $\Delta Z^{\mathrm{QED}}_{12}$ and $\eta_{12}$ are separately gauge independent.

When inserted into the expression for amplitude for the decay $P\to\ell\nu$, the term of order $\alpha_{\mathrm{em}}$ of the renormalised operator $O^{\textrm{W-reg}}_{1}(M_{W})$ of Eq. (81), namely $\delta O^{\textrm{W-reg}}_{1}(M_{W})=O^{\textrm{W-reg}}_{1}(M_{W})-O_{1}^{\chi}$, provides the contribution denoted as $\delta A_{P}^{W}$ in Eq. (30)

[TABLE]

where $X_{P}^{\ell,0}$ is the leptonic trace defined in Eq. (61). We then note that $O_{1}^{\chi}$ and $O_{2}^{\chi}$ entering in Eq. (81) give opposite contributions to the tree-level amplitude, i.e.

[TABLE]

with $A_{P}^{(0)}$ given in Eq. (52) . Therefore, after averaging the amplitude over the values of the parameter $\overline{r}=\pm 1$, in order to cancel out the contribution of the mixing with $O_{3}$ and $O_{4}$, one obtains

[TABLE]

with

[TABLE]

As already noted, the contribution $\delta A_{P}^{W}$ of the matching factor at order $\alpha_{\mathrm{em}}$ to the decay amplitude, expressed by Eqs. (85) and (86), is gauge independent. It then follows that also the order $\alpha_{\mathrm{em}}$ contribution of the bare diagrams to the amplitude, expressed by the other terms in Eq. (30), is by itself gauge independent. Therefore, we can numerically evaluate the two contributions separately by making different choices for the gluon and the photon gauge in the two cases 555It should be noted, however, that while $Z^{\textrm{W-reg}}$ of Eq. (86) is gauge independent at any order of perturbation theory, its actual numerical value may display a residual gauge dependence due to higher order terms in the non-perturbative determination of $\eta_{11}$ which are neglected in the perturbatively evaluated matching coefficient.. In particular, we have chosen to compute the matching factor $Z^{\textrm{W-reg}}$ of Eq. (86) in the Landau gauge for both gluons and photons, because this makes RI' equivalent to RI up to higher orders in the perturbative expansions. On the other hand, in the calculation of the physical amplitudes described in Sec. III (and already computed in Ref. Giusti:2017dwk ) we have used a stochastic photon generated in the Feynman gauge, which has been adopted also in the calculation of $\Gamma^{{\rm pt}}_{0}(L)$ in Ref. Lubicz:2016xro .

As discussed in Ref. Carrasco:2015xwa , when we compute the difference $\Gamma_{0}(L)-\Gamma^{{\rm pt}}_{0}(L)$ in Eq. (1) at leading order in $\alpha_{\mathrm{em}}$, the contribution from the lepton wave function RC cancels out provided, of course, it is evaluated in $\Gamma_{0}(L)$ and $\Gamma^{{\rm pt}}_{0}(L)$ in the same W-regularization scheme and in the same photon gauge. Since $\Gamma^{{\rm pt}}_{0}(L)$ has been computed in Ref. Lubicz:2016xro by omitting the lepton wave function RC contribution in the Feynman gauge, we have to subtract the analogous contribution from Eq. (86) in the Feynman gauge. The QCD and QED corrections to the the lepton wave function RC at $O(\alpha_{\mathrm{em}})$ factorise, so that their contribution does not enter into the non-perturbative determination of the matrix $\eta$, which only contains, by its definition, non-factorisable QCD+QED contributions. Therefore, as discussed in Ref. Carrasco:2015xwa , the subtraction of the lepton wave function RC only requires the replacement of $Z^{\textrm{W-reg}}$ in Eq. (86) by the subtracted matching factor

[TABLE]

where

[TABLE]

The final expression to be used in Eq. (30) is therefore

[TABLE]

with

[TABLE]

To make contact with the factorisation approximation introduced in Refs. Giusti:2017dmp ; Giusti:2017jof , we rewrite Eq. (90) as

[TABLE]

where $Z^{\textrm{W-reg}}_{\eta=0}$ is the result in the factorisation approximation (i.e. with $\eta=0$)

[TABLE]

and $Z^{\mathrm{\,fact}}$ is the factor correcting the result for $\tilde{Z}^{\textrm{W-reg}}$ to include the entries of the matrix $\eta$ determined in Ref. dicarlo

[TABLE]

The values of the coefficients $Z^{\textrm{W-reg}}_{\eta=0}$ and $Z^{\mathrm{fact}}$ are collected in Table 1 for the three values of the inverse coupling $\beta$ adopted in this work and for $\mu=1/a$. In the same Table we also include the values of the coefficient $Z_{m}^{\mathrm{fact}}$ corresponding to the non-factorisable e.m. corrections to the mass RC (see Eq. (40)), evaluated in Ref. dicarlo . The two methods M1 and M2 correspond to different treatments of the $O(a^{2}\mu^{2})$ discretisation effects and are described in Ref. Carrasco:2014cwa . The difference of the results obtained with these two methods enters into the systematic uncertainty labelled as $()_{input}$ in Sec. VI below. The results in Table 1 show that the non-factorisable corrections are significant, of O(12 - 25%) for $Z^{\textrm{W-reg}}$ and even larger, O(40 - 60%), for ${Z}_{m}$.

We close this section by noting that Eq. (29) implies that the contribution to $\delta R_{P}$ from the matching factor in Eq. (89) is $2\tilde{Z}^{\textrm{W-reg}}$. Such a term is mass independent. Thus, as already pointed out in Ref. Giusti:2017dwk , all the matching and mixing contributions to the axial amplitude in Eq. (30) cancel exactly in the difference between the corrections corresponding to two different channels, e.g. in $\delta R_{K}-\delta R_{\pi}$. A similar cancelation also occurs in the difference between the corrections to the amplitudes corresponding to the meson P decaying into two different final-state leptonic channels.

V Finite volume effects at order ${\cal{O}}(\alpha_{\mathrm{em}})$

The subtraction $\Gamma_{0}(L)-\Gamma_{0}^{\textrm{pt}}(L)$ in Eq. (1) cancels both the IR divergences and the structure-independent FVEs, i.e. those of order ${\cal{O}}(1/L)$. The point-like decay rate $\Gamma_{0}^{\textrm{pt}}(L)$ is given by

[TABLE]

where

[TABLE]

with the coefficients $b_{j}$ ($j=\mathrm{IR},0,1,2,3$) depending on the dimensionless ratio $m_{\ell}/M_{P}$ and given explicitly in Eq. (98) of Ref. Lubicz:2016xro (see also Ref. Tantalo:2016vxk ) after the subtraction of the lepton self-energy contribution in the Feynman gauge. An important result of Ref. Lubicz:2016xro is that the structure-dependent FVEs start at order ${\cal{O}}(1/(M_{P}L)^{2})$. Consequently the coefficients $b_{\mathrm{IR},0,1}$ in the factor $Y_{P}^{\ell}(L)$ are “universal”, i.e. they are the same as in the full theory when the structure of the meson $P$ is considered666Notice that the decay rate in the full theory, $\Gamma_{0}(L)$, can be affected also by non-universal FVEs of order ${\cal{O}}[1/(M_{P}L)^{n}]$ with $n\geq 4$ that do not appear in $\Gamma_{0}^{pt}(L)$..

Eq. (30) is therefore replaced by

[TABLE]

where $\delta A_{P}^{W}$ is given by Eq. (89).

In order to study the FVEs in detail we consider four ensembles generated at the same values of $\beta$ and quark masses, but differing in the size of the lattice; these are the ensembles A40.40, A40.32, A40.24 and A40.20 (see Appendix A). The residual FVEs after the subtraction of the universal terms as in Eq. (96) are illustrated in the plots in Fig. 9 for $\delta R_{\pi}$ and $\delta R_{K}$ in the fully inclusive case, i.e. where the energy of the final-state photon is integrated over the full phase space. In this case $\Delta E_{\gamma}=\Delta E_{\gamma}^{\mathrm{max},P}=M_{P}(1-m_{\mu}^{2}/M_{P}^{2})/2$, which corresponds to $\Delta E_{\gamma}^{\mathrm{max},K}\simeq 235$ MeV and $\Delta E_{\gamma}^{\mathrm{max},\pi}\simeq 29$ MeV, respectively. With a muon as the final state lepton, the contribution from photons with energy greater than about 20 MeV is negligible and hence the point-like approximation is valid. In the top plot the universal FV corrections have been subtracted and so we would expect the remaining effects to be of order ${\cal{O}}(1/(M_{P}L)^{2})$ and this is indeed what we see.

In the bottom plot of Fig. 9, in addition to subtracting the universal FVEs, we also subtract the contribution to the order ${\cal{O}}(1/(M_{P}L)^{2})$ corrections from the point-like contribution to $b_{2}$, which can be found in Eq. (3.2) of Ref. Tantalo:2016vxk . We observe that this additional subtraction does not reduce the ${\cal{O}}(1/(M_{P}L)^{2})$ effects, underlining the expectation that these effects are indeed structure dependent.

It can be seen that after subtraction of the universal terms the residual structure-dependent FVEs are almost linear in $1/L^{2}$, which implies that the FVEs of order ${\cal{O}}(1/(M_{P}L)^{3})$ are quite small; indeed they are too small to be resolved with the present statistics. Nevertheless, since the QEDL formulation of QED on a finite box, which is adopted in this work, violates locality Patella:2017fgk , we may expect that there are also FVEs of order ${\cal{O}}(a^{3}/L^{3})$ Tantalo:2016vxk . We have checked explicitly that the addition of such a term in fitting the results shown in Fig. 9 changes the extrapolated value at infinite volume well within the statistical errors.

A more detailed description of the full analysis, including the continuum and chiral extrapolations, is given in the following section. As far as the FVEs are concerned, the central value is obtained by subtracting the universal terms and fitting the residual ${\cal{O}}(1/L^{2})$ corrections to

[TABLE]

where $K_{P}$ and $K_{P}^{\ell}$ are constant fitting parameters and $E_{P}^{\ell}$ is the energy of the charged lepton in the rest frame of the pseudoscalar $P$ (see Eq. (98) below). Such an ansatz is introduced to model the unknown dependance of $b_{2}$ on the ratio $m_{\ell}/M_{P}$. For the four points in each of the plots of Fig. 9 $m_{\ell}/M_{P}$ takes the same value, but this is not true for all the ensembles used in the analysis. We estimate the uncertainty due to the use of the ansatz in Eq. (97) by repeating the same analysis, but on the data in which, in addition to subtracting the universal terms in Eq. (95), we also subtract the term $b_{2}^{\mathrm{pt}}/(M_{P}L)^{2}$, where $b_{2}^{\mathrm{pt}}$ is contribution to $b_{2}$ from a point-like meson Tantalo:2016vxk . Since $b_{2}^{\mathrm{pt}}$ depends on $m_{\ell}/M_{P}$, the result obtained with this additional subtraction is a little different from that obtained with only the universal terms removed and we take the difference as an estimate of the residual FV uncertainty.

VI Results for charged pion and kaon decays into muons

We now insert the various ingredients described in the previous sections into the master formula in Eq. (29) for the decays $\pi^{+}\to\mu^{+}\nu[\gamma]$ and $K^{+}\to\mu^{+}\nu[\gamma]$.

The results for the corrections $\delta R_{\pi}$ and $\delta R_{K}$ are shown in Fig. 10, where the “universal” FSEs up to order ${\cal{O}}(1/L)$ have been subtracted from the lattice data (see the empty symbols) and all photon energies (i.e. $\Delta E_{\gamma}=\Delta E_{\gamma}^{\mathrm{max},P}=M_{P}(1-m_{\mu}^{2}/M_{P}^{2})/2$) are included, since the experimental data on $\pi_{\ell 2}$ and $K_{\ell 2}$ decays are fully inclusive. As already pointed out in section I, structure dependent contributions to real photon emission should be included. According to the ChPT predictions of Ref. Cirigliano:2007ga , however, these contributions are negligible in for both kaon and pion decays into muons, while the same does not hold as well for decays into final-state electrons (see Ref. Carrasco:2015xwa ). This important conclusion needs to be explicitly validated by an ongoing dedicated lattice study of the real photon emission amplitudes in light and heavy P-meson leptonic decays.

The combined chiral, continuum and infinite-volume extrapolations are performed using the following SU(2)-inspired fitting function

[TABLE]

where $m_{ud}=\mu_{ud}/Z_{P}$ and $\mu_{ud}$ is the bare (twisted) mass (see Table 2 in Appendix A below), $E_{P}^{\ell}$ is the lepton energy in the P-meson rest frame, $R_{P}^{(0),(1),(2)}$, $D_{P}$, $K_{P}$ and $K_{P}^{\ell}$ are free parameters. In Eq. (98) the chiral coefficient $R_{P}^{(\chi)}$ is known for both pion and kaon decays from Ref. Bijnens:2006mk ; in QED the coefficients are

[TABLE]

while in qQED they are

[TABLE]

where $X$ is obtained from the chiral limit of the ${\cal{O}}(\alpha_{\mathrm{em}})$ correction to $M_{\pi^{\pm}}^{2}$ (i.e. $\delta M_{\pi^{\pm}}^{2}=4\pi\alpha_{\mathrm{em}}Xf_{0}^{2}+{\cal{O}}(m_{ud})$). In Ref. Giusti:2017dmp we found $X=0.658(40)$.

Using Eq. (98) we have fitted the data for $\delta R_{\pi}$ and $\delta R_{K}$ using a $\chi^{2}$-minimization procedure with an uncorrelated $\chi^{2}$, obtaining values of $\chi^{2}/\mbox{d.o.f.}$ always around $0.9$. The uncertainties on the fitting parameters do not depend on the $\chi^{2}$-value, because they are obtained using the bootstrap samplings of Ref. Carrasco:2014cwa (see Appendix A). 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. 10. It can be seen that the residual SD FVEs are still visible in the data and well reproduced by our fitting ansatz in Eq. (98). Discretisation effects on the other hand, only play a minor role.

At the physical pion mass in the continuum and infinite-volume limits we obtain

[TABLE]

where

• •

$()_{stat+fit}$ indicates the uncertainty induced by the statistical Monte Carlo errors of the simulations and its propagation in the fitting procedure;

• •

$()_{input}$ is the error coming from the uncertainties of the input parameters of the quark-mass analysis of Ref. Carrasco:2014cwa ;

• •

$()_{chiral}$ is the difference between including or excluding the chiral logarithm in Eq. (98), i.e. taking $R_{\chi}\neq 0$ or $R_{\chi}=0$;

• •

$()_{FVE}$ is the difference between the analyses of the data corresponding to the FVE subtractions up to the order ${\cal{O}}(1/L)$ alone or by also subtracting the term proportional to $b_{2}^{\mathrm{pt}}/(M_{P}L)^{2}$ (see Fig. 9 and the discussion towards the end of Sec. V);

• •

$()_{disc}$ is the uncertainty coming from including ($D\neq 0$) or excluding (setting $D=0$) the discretisation term proportional to $a^{2}$ in Eq. (98);

• •

$()_{qQED}$ is our estimate of the uncertainty of the QED quenching. This is obtained using the ansatz (98) with the coefficient $R_{\chi}$ of the chiral log fixed either at the value (100), which corresponds to the qQED approximation, or at the value (99), which includes the effects of the up, down and strange sea-quark charges Bijnens:2006mk . The change both in $\delta R_{\pi}^{\textrm{phys}}$ and in $\delta R_{K}^{\textrm{phys}}$ is $\simeq 0.0003$, which has been already added in the central values given by Eqs. (101) and (102). To be conservative, we use twice this value for our estimate of the qQED uncertainty.

Our results in Eqs. (101) - (102) can be compared with the ChPT predictions $\delta R_{\pi}^{\textrm{phys}}=0.0176(21)$ and $\delta R_{K}^{\textrm{phys}}=0.0064(24)$ obtained in Ref. Cirigliano:2011tm and adopted by the PDG  PDG ; Rosner:2015wva . The difference is within one standard deviation for $\delta R_{\pi}^{\textrm{phys}}$, while it is larger for $\delta R_{K}^{\textrm{phys}}$. Note that the precision of our determination of $\delta R_{\pi}^{\textrm{phys}}$ is comparable to the one obtained in ChPT, while our determination of $\delta R_{K}^{\textrm{phys}}$ has a much better accuracy compared to that obtained using ChPT; the improvement in precision is a factor of about 2.2. We stress that the level of precision of our pion and kaon results depends crucially on the non-perturbative determination of the chirality mixing, carried out in section IV by including simultaneously QED at first order and QCD at all orders.

As already stressed, the correction $\delta R_{P}$ and the QCD quantity $f_{P}^{(0)}$ separately depend on the prescription used for the separation between QED and QCD corrections Gasser:2010wz . Only the product $f_{P}^{(0)}\sqrt{1+\delta R_{P}}$ is independent of the prescription and its value, multiplied by the relevant CKM matrix element, yields the P-meson decay rate. We remind the reader that our results (101) - (102) are given in the GRS prescription (see the dedicated discussion in sections II.2 and III) in which the renormalised couplings and quark masses in the full theory and in isosymmetric QCD coincide in the $\overline{\rm MS}$ scheme at a scale of 2 GeV Gasser:2003hk . We remind the reader that, to the current level of precision, this GRS scheme can be considered equivalent to the FLAG scheme.

Taking the experimental values $\Gamma(\pi^{-}\to\mu^{-}\bar{\nu}_{\mu}[\gamma])=3.8408(7)\cdot 10^{7}$ s*-1* and $\Gamma(K^{-}\to\mu^{-}\bar{\nu}_{\mu}[\gamma])=5.134(11)\cdot 10^{7}$ s*-1* from the PDG PDG and using our results (101) - (102), we obtain

[TABLE]

where the first error is the experimental uncertainty and the second is that from our theoretical calculations. The result for the pion in Eq. (103) agrees within the errors with the updated value $f_{\pi}^{(0)}|V_{ud}|=127.12(13)$ MeV PDG , obtained by the PDG and based on the model-dependent ChPT estimate of the e.m. corrections from Ref. Cirigliano:2011tm . Our result for the kaon in Eq. (104) however, is larger than the corresponding PDG value $f_{K}^{(0)}|V_{us}|=35.09(5)$ MeV PDG , based on the ChPT calculation of Ref. Cirigliano:2011tm , by about $2$ standard deviations.

As anticipated in the Introduction and discussed in detail in Sec. III, we cannot use the result (103) to determine the CKM matrix element $|V_{ud}|$, since the pion decay constant was used by ETMC Carrasco:2014cwa to set the lattice scale in isosymmetric QCD and its value, $f_{\pi}^{(0)}=130.41(20)$ MeV, was based on the determination of $|V_{ud}|$ obtained from super-allowed $\beta$-decays in Ref. Hardy:2014qxa . On the other hand, adopting the best lattice determination of the QCD kaon decay constant, $f_{K}^{(0)}=156.11(21)$ MeV FLAG ; Dowdall:2013rya ; Carrasco:2014poa ; Bazavov:2017lyh  777The average value of $f_{K^{\pm}}$ quoted by FLAG FLAG includes the strong IB corrections. In order to obtain $f_{K}^{(0)}$ therefore, we have subtracted this correction which is given explicitly in Refs. Dowdall:2013rya ; Carrasco:2014poa ; Bazavov:2017lyh ., we find that Eq. (104) implies

[TABLE]

which is a result with the excellent precision of $\simeq 0.2\%$.

Since the non-factorisable e.m. corrections to the mass RC (see the coefficient $Z_{m}^{\mathrm{fact}}$ in Table 1) were not included in Ref. Giusti:2017dwk , we update our estimate of the ratio of the kaon and pion decay rates:

[TABLE]

Using the pion and kaon experimental decay rates we get

[TABLE]

Using the best $N_{f}=2+1+1$ lattice determination of the ratio of the QCD kaon and pion decay constants, $f_{K}^{(0)}/f_{\pi}^{(0)}=1.1966~{}(18)$ FLAG ; Dowdall:2013rya ; Carrasco:2014poa ; Bazavov:2017lyh , we find

[TABLE]

Taking the updated value $|V_{ud}|=0.97420\,(21)$ from super-allowed nuclear beta decays Hardy:2016vhg , Eq. (108) yields the following value for the CKM element $|V_{us}|$:

[TABLE]

which agrees with our result (105) within the errors. Note that our result (109) agrees with the latest estimate $|V_{us}|=0.2253(7)$, recently updated by the PDG PDG , but it improves the error by a factor of approximately 1.5.

Taking the values $|V_{ub}|=0.00413(49)$ PDG and $|V_{ud}|=0.97420(21)$ Hardy:2016vhg our result in Eq. (109) implies that the unitarity of the first-row of the CKM matrix is confirmed to better than the per-mille level

[TABLE]

With the same value $|V_{ud}|=0.97420(21)$ from super-allowed nuclear beta decays Hardy:2016vhg , our result (103) implies for the QCD pion decay constant (in the GRS prescription) the following value

[TABLE]

which, as anticipated in Sec. III, agrees within the errors with the value $f_{\pi}^{(0)}=130.41\,(20)$ MeV adopted in Ref. Carrasco:2014cwa to set the lattice scale in the isosymmetric QCD theory. This demonstrates the equivalence of the GRS and PDG schemes within the precision of our simulation.

In a recent paper Seng:2018yzq the hadronic contribution to the electroweak radiative corrections to neutron and super-allowed nuclear $\beta$ decays has been analyzed in terms of dispersion relations and neutrino scattering data. With respect to the result $V_{ud}=0.97420(21)$ from Ref. Hardy:2016vhg a significant shift in the central value and a reduction of the uncertainty have been obtained, namely $V_{ud}=0.97370(14)$ Seng:2018yzq . The impact of the new value of $V_{ud}$ on our determinations of $V_{us}$ and $f_{\pi}^{(0)}$ is $V_{us}=0.22526\,(46)$ and $f_{\pi}^{(0)}=130.72\,(12)$ MeV, i.e. well within the uncertainties shown in Eqs. (109) and (111), respectively. On the contrary, the first-row CKM unitarity (110) will be significantly modified into

[TABLE]

which would imply a $\simeq 3.4\sigma$ tension with unitarity. A confirmation of the new calculation of the radiative corrections made in Ref. Seng:2018yzq is therefore urgently called for.

Before closing this section, we comment briefly about the comparison between our result $\delta R_{K}^{\mathrm{phys}}=0.0024(10)$ and the corresponding model-dependent ChPT prediction $\delta R_{K}^{\mathrm{phys}}=0.0064(24)$ from Ref. Cirigliano:2011tm . The latter is obtained by adding a model-dependent QED correction of $0.0107(21)$ and a model-independent next-to-leading strong IB contribution equal to $-0.0043(12)$. Our result on the other hand, obtained in the GRS prescription, stems from a QED correction equal to $0.0088(9)$ and a strong IB term equal to $-0.0064(7)$ (see also Ref. Giusti:2017xrv ). The difference between our result and the ChPT prediction of Ref. Cirigliano:2011tm appears to be mainly due to a different strong IB contribution. Thus, in the present $N_{f}=2+1+1$ study, we confirm for the strong IB term a discrepancy at the level of about 2 standard deviations, which was already observed at $N_{f}=2$ in Ref. deDivitiis:2013xla .

VII Conclusions

In this paper we have presented the details of the first lattice computation of the leading e.m. and strong IB corrections to the $\pi^{+}\to\mu^{+}\nu$ and $K^{+}\to\mu^{+}\nu$ leptonic decay rates, following a method recently proposed in Ref. Carrasco:2015xwa . This expands significantly on the discussion of Ref. Giusti:2017dwk , where the results and a brief outline of the calculation had been presented. The results were 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.

The effective weak Hamiltonian in the W-regularisation scheme appropriate for this calculation is obtained from the bare lattice operators in two stages. First of all, the lattice operators are renormalised non-perturbatively in the RI′-MOM scheme at $O(\alpha_{\mathrm{em}})$ and to all orders in the strong coupling $\alpha_{s}$. Because of the breaking of chiral symmetry in the twisted mass formulation we have adopted this renormalisation includes the mixing with other four-fermion operators of different chirality. In the second step we perform the matching from the RI′-MOM scheme to the W-regularisation scheme perturbatively. By calculating and including the two-loop anomalous dimension at $O(\alpha_{\mathrm{em}}\alpha_{s})$ Brod:2008ss , the residual truncation error of this matching is of $O(\alpha_{\mathrm{em}}\alpha_{s}(M_{W}))$, reduced from $O(\alpha_{\mathrm{em}}\alpha_{s}(1/a))$ in our earlier work Giusti:2017dwk ; Carrasco:2015xwa .

The evaluation of isospin breaking (IB) “corrections” raises the question of how QCD without these corrections is defined. Since IB corrections change hadronic masses and other physical quantities a prescription is needed to define QCD, whether isosymmetric or not, and in Sec. II and Appendix B we discuss this issue in detail. In particular, the correction $\delta R_{P}$ and the QCD quantity $f_{P}^{(0)}$ separately depend on the prescription used for the definition of QCD Gasser:2010wz . Only the product $f_{P}^{(0)}\sqrt{1+\delta R_{P}}$ is independent of the prescription and its value, multiplied by the relevant CKM matrix element, yields the P-meson decay rate. In this paper we chose to follow the conventionally used GRS prescription (see the dedicated discussion in sections II.2 and III) in which the renormalised couplings and quark masses in the full QCD+QED theory and in isosymmetric QCD coincide in the $\overline{\rm MS}$ scheme at a scale of 2 GeV Gasser:2003hk . For future studies however, we advocate the use of “hadronic schemes” in which QCD is defined by requiring that a set of hadronic quantities (for example a set of hadronic masses) take their physical values in QCD and in QCD+QED.

The main results of the calculation are presented in Sec. VI together with a detailed discussion of their implications. In summary, after extrapolation of the data to the physical pion mass, and to the continuum and infinite-volume limits, the isospin-breaking corrections to the leptonic decay rates can be written in the form:

[TABLE]

where $\Gamma^{(0)}$ is the leptonic decay rate at tree level in the GRS scheme (see Eqs. (101) and (102)). These results can be compared with the ChPT predictions $\delta R_{\pi}^{\textrm{phys}}=0.0176(21)$ and $\delta R_{K}^{\textrm{phys}}=0.0064(24)$ obtained in Ref. Cirigliano:2011tm and adopted by the PDG  PDG ; Rosner:2015wva . The difference is within one standard deviation for $\delta R_{\pi}^{\textrm{phys}}$, while it is larger for $\delta R_{K}^{\textrm{phys}}$. We also underline that our result $|V_{us}|=0.22538(46)$ in Eq. (109), together with the value of $V_{ud}$ determined in Ref. Hardy:2016vhg and $|V_{ub}|$ from the PDG PDG , implies that the unitarity of the first row of the CKM matrix is satisfied at the per-mille level (see Eq. (110)).

Acknowledgments

We gratefully acknowledge the CPU time provided by the Partnership for Advanced Computing in Europe (PRACE) under the project Pra10-2693 “QED corrections to meson decay rates in Lattice QCD” and by CINECA under the specific initiative INFN-LQCD123 on the BG/Q system Fermi at CINECA (Italy). V.L., G.M. and S.S. thank MIUR (Italy) for partial support under the contract PRIN 2015. G.M. also acknowledges partial support from the European Research Council Ideas Advanced Grant n. 267985 “DaMeSyFla”. C.T.S. was partially supported by the UK Science and Technology Facilities Council (STFC) grant No. ST/P000711/1 and by an Emeritus Fellowship from the Leverhulme Trust. N.T. gratefully acknowledges the University of Rome Tor Vergata for the support granted to the project PLNUGAMMA.

Appendix A Details of the simulation

The gauge ensembles used in this work are those generated by ETMC with $N_{f}=2+1+1$ dynamical quarks and used in Ref. Carrasco:2014cwa to determine the up, down, strange and charm quark masses. We use the Iwasaki action Iwasaki:1985we for the gluons and the Wilson Twisted Mass Action Frezzotti:2000nk ; Frezzotti:2003xj ; Frezzotti:2003ni for the sea quarks. In the valence sector we adopt a non-unitary setup Frezzotti:2004wz in which the strange quark is regularized as an Osterwalder-Seiler fermion Osterwalder:1977pc , while the up and down quarks have the same action as the sea. Working at maximal twist such a setup guarantees an automatic ${\cal{O}}(a)$-improvement Frezzotti:2003ni ; Frezzotti:2004wz .

We have performed simulations at three values of the inverse bare lattice coupling $\beta$ and at several different lattice volumes as shown in Table 2. We allow a separation of 20 trajectories between each of the $N_{\mathrm{cfg}}$ analysed configurations. For the earlier investigation of finite-volume effects (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 such an investigation, which is crucial in the present work, we have generated a further gauge ensemble, A40.40, at a larger value of the lattice size $L$.

At each lattice spacing, different values of the light sea-quark masses have been considered. The light valence and sea quark masses are always taken to be degenerate. The bare mass of the valence strange quark ($a\mu_{s}$) is obtained, at each $\beta$, using the physical strange mass and the mass RCs determined in Ref. Carrasco:2014cwa . There the “FLAG” hadronic scheme was adopted in which the pion and kaon masses in isosymmetric QCD are equal to $M_{\pi}^{(0),\textrm{FLAG}}=134.98$ MeV and $M_{K}^{(0),\textrm{FLAG}}=494.2$ MeV and the lattice scale is fixed by the value $f_{\pi}^{(0),\textrm{FLAG}}=130.41~{}(20)$ MeV for the physical pion decay constant. In the charm sector instead, the $D_{s}$-meson mass $M_{D_{s}}^{(0)}$ was chosen to be equal to its experimental value $M_{D_{s}^{+}}=1969.0(1.4)$ MeV PDG . The values of the lattice spacing are found to be $a=0.0885(36)$, $0.0815(30)$, $0.0619(18)$ fm at $\beta=1.90$, $1.95$ and $2.10$, respectively.

The two valence quarks $q_{1}$ and $q_{2}$ in the P meson are regularized with opposite values of the Wilson $r$-parameter ($r_{2}=-r_{1}$) in order to guarantee that discretisation effects on the P meson mass are of order ${\cal{O}}(a^{2}\mu~{}\Lambda_{\textrm{QCD}})$. The lepton is a free twisted-mass fermion with Wilson parameter $r_{\ell}$ and its mass is taken fixed at the physical muon value $m_{\ell}=m_{\mu}=105.66$ MeV  PDG . The regularization of the (massless) neutrino is irrelevant and it is taken to be a free fermion field.

In this work we made use of the bootstrap samples generated for the input parameters of the quark mass analysis of Ref. Carrasco:2014cwa . There, eight branches of the analysis were adopted differing in:

• •

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;

• •

the chiral extrapolation performed with fitting functions chosen to be either a polynomial expansion or a Chiral Perturbation Theory (ChPT) Ansatz in the light-quark mass;

• •

the choice between the methods M1 and M2, which differ by ${\cal{O}}(a^{2})$ effects, used to determine the mass RC $Z_{m}=1/Z_{P}$ in the RI′-MOM scheme .

Appendix B Relating observables in the full theory and in QCD

In this appendix we provide the detailed derivation of the relation between observables calculated in the full theory (QCD+QED) and in QCD (in the absence of QED). We start by a discussion of the separation of the QCD action from that in the full theory.

B.1 Actions of the full theory and of QCD

The lattice action in the full theory given in Eq. (4) can be written as

[TABLE]

where $S_{\ell,0}=S^{\textrm{kin}}_{\ell,0}+m_{\ell}S_{\ell}^{m}$ and the counterterm $S^{\textrm{ct}}$ and the QED vertices $\Delta S$ are given by

[TABLE]

We now consider these terms in detail using Wilson fermions for illustration. The kinetic term for the quark with flavour $f$, $S_{f}^{\textrm{kin}}$, is given by:

[TABLE]

where $\psi_{f}$ is the quark field, while $U_{\mu}$ and $V_{f,\mu}$ are the QCD and QED gauge links, respectively. Specifically

[TABLE]

where $e_{f}$ is the charge of the quark with flavour $f$ in units of the positron charge. The forward and backward derivatives are given by

[TABLE]

The leptonic action is given by

[TABLE]

with $\psi_{\ell}$ being the lepton field. The renormalisation of the lepton masses is performed perturbatively, by requiring that the on-shell masses correspond to the physical ones.

In QCD the kinetic term only includes the gluon links so that for Wilson fermions

[TABLE]

and the derivatives are defined in Eqs. (120) - (121) with the $V_{f}=1$. Since for leptonic and semileptonic decays, leptonic spinors are present even in the absence of electromagnetism, it is also convenient to define the kinetic action for free leptons:

[TABLE]

B.2 Relation between observables in the full theory and in QCD

Physical observables are determined from correlation functions evaluated from lattice computations in the full theory. For a generic observable ${\cal{O}}$ evaluated in the full theory up to ${\cal{O}}(\alpha_{\mathrm{em}})$ we write:

[TABLE]

where in the integrand $O$ is a multilocal composite operator. For a given choice of the strong coupling $g_{s}$, the parameters of the action, the bare quark masses and the lattice spacing, are determined by imposing that a set of physical quantities take their experimental values as explained in Sec. II.1. Physical quantities other than those used for the calibration can now be determined unambiguously up to lattice artefacts, which are removed by taking the continuum limit.

In general, the determination of physical observables requires the processing of correlation functions of the form of Eq. (125). Hadronic masses, for example are obtained from the behaviour in the time separation of two interpolating operators and the determination of hadronic matrix elements may require the cancelation of interpolating operators at the source and/or sink using a combination of three- and two-point functions. The discussion in this appendix concerns the evaluation of a generic correlation function.

We now turn to the definition of correlation functions in QCD defined in a generic scheme. For a generic observable ${\cal{O}}$ we define its value in QCD by:

[TABLE]

where the bare quark masses and the lattice spacing are defined as discussed in Sec. II.2. The free QED action is included in the numerator and denominator of Eq. (126) since even without radiative corrections the physical quantities such as $\Gamma(K_{\ell 2})$ and $\Gamma(\pi_{\ell 2})$ studied in this paper are obtained by combining the results for hadronic matrix elements obtained from QCD simulations with leptonic spinors. Moreover, for other quantities, for example the long-distance contributions to the amplitude for the rare kaon decay $K^{+}\to\pi^{+}\nu\bar{\nu}$, there are internal free lepton propagators even in the absence of isospin breaking Christ:2016eae .

Comparing Eqs. (125) and (126) we arrive at

[TABLE]

where the subscript “conn” reminds that only connected Feynman diagrams contribute: $\langle O_{1}O_{2}\rangle_{\textrm{conn}}=\langle O_{1}O_{2}\rangle-\langle O_{1}\rangle\,\langle O_{2}\rangle$.

There is one final subtlety which we must account for. We need to convert the results obtained from simulations in lattice units (i.e. in units of the lattice spacing) into values given in physical units such as MeV. Eq. (127) is also written in lattice units. Imagine that the observable $O$ has mass dimension $n$ and rewrite Eq. (127) with the lattice spacing included explicitly:

[TABLE]

where, since we are working to first order in isospin breaking, in the second term on the right-hand side we do not need to distinguish between the lattice spacing in the full theory ($a$) and that obtained in QCD ($a_{0}$). The quantity which we wish to determine, $\langle O\rangle^{\textrm{full}}$ in physical units, is therefore given by

[TABLE]

where $\delta a=a-a_{0}$. The three expectation values on the right-hand side of (129) are directly computed in QCD simulations.

Appendix C Non-perturbative renormalisation in the RI'-MOM scheme

In this paper, as explained in Sec. IV, we have renormalised the weak four-fermion operator $O_{1}$ non-perturbatively on the lattice to all orders in $\alpha_{s}$ and up to first order in $\alpha_{\mathrm{em}}$. In this appendix we describe the main steps of the non-perturbative renormalisation procedure at ${\cal O}(\alpha_{\mathrm{em}})$ and we refer the reader to a forthcoming publication dicarlo for further details and results.

Given the amputated Green function, $\Lambda_{O}$, of an operator $O$ computed in a given gauge between external states with momentum $p$ and a suitable projector on the relevant Dirac structure, $P_{O}$, we define the projected Green function as

[TABLE]

In the RI'-MOM scheme, the renormalisation constant (RC) $Z_{O}(\mu a)$ is found by imposing the condition Martinelli:1994ty

[TABLE]

where

[TABLE]

The $Z_{f}$ are the RCs of the external fields and the index $f$ runs over all external fields entering the expression of the composite operator $O$. For the four-fermion operators considered in this work, the RCs $Z_{O}(\mu a)$ and the projected Green functions $\Gamma_{O}(pa)$ are $5\times 5$ matrices, the latter with elements $(\Gamma_{O})_{ij}=\mathrm{Tr}\left[\Lambda_{O_{i}}\,P_{O_{j}}\right]$. In QCD+QED the RCs $Z_{O}$ and $Z_{f}$ depend both on the strong and the e.m. coupling constants.

Following the discussion of Sec. IV (see Eqs. (75) - (78)), we write the RCs of any composite operator, and in particular of the fields, bilinear and four-fermion operators, in the generic decomposition

[TABLE]

where $Z_{O}^{\mathrm{QCD}}$ and $Z_{O}^{\mathrm{QED}}$ are the RCs of the operator $O$ in pure QCD and pure QED respectively and we have put

[TABLE]

The first term, $\Delta Z_{O}^{\mathrm{QED}}$, in Eq. (134) represents the pure QED contribution to the RC at ${\cal O}(\alpha_{\mathrm{em}})$, whereas $\eta_{O}$ contains the ${\cal O}(\alpha_{\mathrm{em}})$ non-factorisable QCD+QED correction.

In terms of the QCD renormalised operators $O^{\chi}$, as those in Eq. (79), we define the QCD renormalised projected Green function $\Gamma_{O}^{\chi}$, and expand it at first order in $\alpha_{\mathrm{em}}$

[TABLE]

where we have used the RI'-MOM renormalisation condition $Z_{\Gamma_{O}}^{\mathrm{QCD}}(\mu a)\Gamma^{\mathrm{QCD}}_{O}(\mu a)=1$ applied in the pure QCD theory and defined

[TABLE]

Using Eqs. (133) and (135), we can rewrite Eq. (131) at first order in $\alpha_{\mathrm{em}}$ as

[TABLE]

which provides, in turn, the RI-MOM renormalisation condition at order $\alpha_{\mathrm{em}}$

[TABLE]

Using the expression of $Z_{\Gamma_{O}}$ in Eq. (132) in terms of $Z_{O}$ and the external fields RCs one also obtains

[TABLE]

Thus, $\Delta Z_{O}$ is expressed directly in terms of the ${\cal O}(\alpha_{\mathrm{em}})$ contribution to the QCD renormalised projected Green function $\Delta\Gamma_{O}^{\chi}=Z_{\Gamma_{O}}^{\mathrm{QCD}}\Delta\Gamma_{O}$ evaluated at $p^{2}=\mu^{2}$.

In the following we describe a completely non-perturbative determination of the RCs $\Delta Z_{O}(\mu a)$ to all orders in $\alpha_{s}$. We will assume that all the relevant RCs of fields and composite operators in pure QCD have been already determined, by following the standard RI'-MOM renormalisation procedure. With appropriate modifications to the kinematical conditions and projectors, the discussion can readily be adapted to similar schemes, such as the Symmetric Momentum subtraction one Sturm:2009kb .

In addition to the renormalisation of the four-fermion operator appearing in the Hamiltonian, the e.m. shift of the quark masses (see Sec. III.1) requires the knowledge of the RC of the pseudoscalar density deDivitiis:2013xla . We therefore start by discussing the non-perturbative renormalisation of quark bilinear operators.

C.1 Renormalisation of the quark field and bilinear operators

We start with the renormalisation of the quark fields. The e.m. corrections to a quark propagator can be represented schematically in the form

[TABLE]

where the last two diagrams represent the mass and critical Wilson parameter counter-terms deDivitiis:2013xla .

The amputated one-particle irreducible two-point function is then given by

[TABLE]

and the correction to the quark field RC in the RI'-MOM scheme is obtained, according to Eq. (138), as

[TABLE]

The e.m. correction to the RC $Z_{O}$ of a generic bilinear operator $O_{\Gamma}=\bar{q}_{2}\,\Gamma\,q_{1}$, where $\Gamma$ is one of the Dirac matrices (${\Gamma=1,\gamma^{5},\gamma^{\mu},\gamma^{\mu}\gamma^{5},\sigma^{\mu\nu}}$), is given by Eq. (139), which in this case reads

[TABLE]

Two kinds of corrections contribute to the amputated Green function: either the QCD Green function is amputated with the e.m. corrections on the inverse propagators, or the correction to the Green function itself is amputated with QCD propagators. Thus we have

[TABLE]

with

[TABLE]

where $G_{O}$ is the non-amputated Green function and $\Delta G_{O}$ is given diagrammatically by

[TABLE]

In this work we have used an improved method to compute the first diagram in Eq. (146), as well as all the diagrams containing a photon propagator connecting different points. In this method, some of the sequential propagators introduced in Ref. Giusti:2017dmp are summed in order to reduce the number of inversions of the Dirac matrix. All details of the calculation will be given in the forthcoming publication dicarlo .

Before closing this subsection, we stress that in the calculation of $Z_{P}$ and its e.m. correction $\Delta Z_{P}$, the Goldstone pole contamination has been taken into account and subtracted. In pure QCD, at each $p^{2}$ and for each combination of valence quark masses, $\mu_{1}$ and $\mu_{2}$, the amputated Green function has been fitted to the ansatz

[TABLE]

where $M_{P}\equiv M_{P}(\mu_{1},\mu_{2})$ is the mass of the pseudoscalar meson composed of valence quarks of mass $\mu_{1}$ and $\mu_{2}$. When including QED in the calculation, Eq. (147) has to be modified to take into account the e.m. correction to the meson mass. By considering the ansatz in Eq. (147) in QCD+QED and expanding it in terms of $\alpha_{\mathrm{em}}$ one finds

[TABLE]

where $\Delta M_{P}^{2}$ is the correction to $M_{P}^{2}$ evaluated in Ref. Giusti:2017dmp . Note, in particular, that $\Delta\Gamma_{P}$ also receives the contribution of a double pole. In Eq. (148) only the coefficients $A_{1}$, $B_{1}$, $C_{1}$ need to be fitted, since the values of $B_{0}$ and $C_{0}$ are already obtained from the QCD fit in Eq. (147).

C.2 Renormalisation of the four-fermions operators

We conclude this section by describing the calculation of the RCs of the complete basis of four-fermion operators $O_{i}$ ($i=1,\ldots,5$), in the RI'-MOM scheme. In this case, the renormalisation condition (139) for the renormalisation matrix at ${\cal{O}}(\alpha_{\mathrm{em}})$ reads:

[TABLE]

where $\Delta Z_{\ell}$ is only e.m. and can be computed in perturbation theory. We remind the reader that this term is omitted in the actual calculation since its contribution cancels out in the difference $\Gamma_{0}(L)-\Gamma_{0}^{\textrm{pt}}(L)$.

In Eq. (149) $\Delta\Gamma_{O}^{\chi}$ is a matrix expressed by

[TABLE]

As in the case of bilinear operators, the correction to the amputated Green function gets two kind of contributions,

[TABLE]

and in this case ${\Delta G}_{O_{i}}$ is given by

[TABLE]

The fermionic lines on the left-hand side of the diagrams in Eq. (C.2) represent the ingoing and outgoing light quarks. On the right-hand side, the external charged anti-lepton and the neutrino propagators are drawn for illustration but not actually included in the calculation. For this reason their amputation is neglected in Eq. (151). The lepton self-energy is not reported in Eq. (C.2) since its contribution cancels out in the amputation.

Appendix D Matching, chirality mixing and fermion operators in the twisted mass regularisation

In the main text and in Appendix C we have described the renormalisation of the relevant operators in the physical basis. This discussion is valid for a generic Wilson-like fermion regularisation. In this appendix we address instead some important aspects peculiar to the twisted mass fermions used in our numerical calculation. We derive, in particular, the relations between RCs in the so-called physical and twisted basis, for the bilinear and four-fermion operators considered in this work.

The relevant observation is that the lattice action for twisted mass fermions at maximal twist in the twisted basis only differs from the standard Wilson fermion lattice action for the twisted rotation of the fermion mass term. The two actions become identical in the chiral limit. It then follows that, in any mass-independent renormalisation scheme, the RCs for twisted mass operators in the twisted basis are the same as those of the corresponding operators with standard Wilson fermions. It is customary to denote these RCs, for a generic operator $O$, as $Z_{O}$. They are valid for both standard Wilson and twisted mass operators in the twisted basis and differ, in general, from the RCs for twisted mass operators in the physical basis, that we denote here as $Z_{O}^{(0)}$.

At maximal twist the rotation from the twisted to the physical basis for both quark and lepton fields is given by

[TABLE]

where $q$ and $\ell$ are the quark and lepton fields in the physical basis and $r_{q}$ and $r_{\ell}$ are the corresponding $r$-parameters. In our simulations we use opposite values of the $r$-parameter for the two valence quarks, ${r_{2}=-r_{1}}$ ($r_{i}=\pm 1$). The quark and lepton bilinears then transform as

[TABLE]

From Eqs. (D) one readily derives the relations between the quark vector and axial vector current in the two basis,

[TABLE]

which, in turn, determine the relation between the RCs in the two basis

[TABLE]

where $\hat{O}$ denotes the generic renormalised operator. One then sees from Eq. (156) that the RC $Z_{V}^{(0)}$ of the vector current in the physical basis, with $r_{1}=-r_{2}$, is simply the RC of the axial current in the twisted basis, which in turn is just $Z_{A}$ computed with Wilson fermions in the chiral limit. Analogously, $Z_{A}^{(0)}$ in the physical basis, with $r_{1}=-r_{2}$, corresponds to $Z_{V}$ computed with Wilson fermions in the chiral limit.

From the transformations (D) one can also derive the relations between the four-fermion operators $O_{1}-O_{5}$ of Eq. (24) and (73) in the physical and twisted basis

[TABLE]

We can then obtain the relation between the renormalisation matrix in the physical basis, $Z^{(0)}$, and the corresponding matrix $Z$ for standard Wilson fermions. In particular, for the weak operator $O_{1}$ one finds

[TABLE]

with $\overline{r}\equiv r_{1}r_{\ell}$. Therefore

[TABLE]

Eq. (159) shows in particular that the mixing coefficients $Z_{13,14,15}$ for the operators $O_{3,4,5}$ are proportional to $\overline{r}\equiv r_{1}r_{\ell}$. Thus, we can eliminate the mixing with these operators by simply averaging the numerical results over the two possible values $\overline{r}=\pm 1$.

In order to illustrate the above point, using the results of Ref. Carrasco:2015xwa obtained in perturbation theory at order ${\cal{O}}(\alpha_{s}^{0})$, the coefficients $\Delta Z^{QED}_{1j}=Z^{(0)}_{1j}/(\alpha_{\mathrm{em}}/4\pi)$ are explicitly given in the physical basis by

[TABLE]