First lattice calculation of radiative leptonic decay rates of pseudoscalar mesons

We present a non-perturbative lattice calculation of the form factors which contribute to the amplitudes for the radiative decays $P\to \ell \bar \nu_\ell \gamma$, where $P$ is a pseudoscalar meson and $\ell$ is a charged lepton. Together with the non-perturbative determination of the corrections to the processes $P\to \ell \bar \nu_\ell$ due to the exchange of a virtual photon, this allows accurate predictions at $O(\alpha_{em})$ to be made for leptonic decay rates for pseudoscalar mesons ranging from the pion to the $D_s$ meson. We are able to separate unambiguously and non-pertubatively the point-like contribution, from the structure-dependent, infrared-safe, terms in the amplitude. The fully non-perturbative $O(a)$ improved calculation of the inclusive leptonic decay rates will lead to the determination of the corresponding Cabibbo-Kobayashi-Maskawa (CKM) matrix elements also at $O(\alpha_{em})$. Prospects for a precise evaluation of leptonic decay rates with emission of a hard photon are also very interesting, especially for the decays of heavy $D$ and $B$ mesons for which currently only model-dependent predictions are available to compare with existing experimental data.


I. INTRODUCTION
The unitarity of the CKM matrix is one of the most precise tests of the Standard Model.
Indeed, CKM unitarity may rule out many theoretically well motivated models for new physics and put severe constraints on the energy scale where new phenomena might occur, well beyond the range accessible to direct experimental searches. In this respect, leptonic decay rates of light and heavy pseudoscalar mesons are essential ingredients for the extraction of the CKM matrix elements. A first-principles calculation of these quantities requires non-perturbative accuracy and hence numerical lattice simulations. Moreover, in order to fully exploit the presently available experimental information and to perform the next generation of flavour-physics tests, O(α em ) electromagnetic corrections must be included. In this endeavour, the radiative leptonic decays P → ν (γ) (where P is a negatively charged pseudoscalar meson, a lepton,ν the corresponding anti-neutrino and γ a photon) are particularly important, see [1].
Knowledge of the radiative leptonic decay rate in the region of small (soft) photon energies is required in order to properly define the infrared-safe measurable decay rate for the process P → ν (γ). Indeed, according to the well-known Bloch-Nordsieck mechanism [2], the integral of the radiative decay rate in the phase space region corresponding to soft photons must be added to the decay rate with no real photons in the final states (the so-called virtual rate) in order to cancel infrared divergent contributions appearing in unphysical quantities at intermediate stages of the calculations.
On the one hand, in the limit of ultra-soft photon energy the radiative decay rate can be reliably calculated in an effective theory in which the meson is treated as a point-like particle. This is a manifestation of the well-known mechanism known as the "universality of infrared divergences" (see for example Ref. [3,4]) that finds its physical explanation in the fact that ultra-soft photons cannot resolve the internal structure of the meson. On the other hand, the ultra-soft limit is an idealisation and experimental measurements, particularly in the case of heavy mesons, are inclusive up to photon energies that may be too large to safely neglect the Structure-Dependent (SD) corrections to the point-like approximation.
In the region of hard (experimentally detectable) photon energies, radiative leptonic decays represent important probes of the internal structure of the mesons. Moreover, radiative decays can provide independent determinations of CKM matrix elements with respect to the purely leptonic channels. A non-perturbative calculation of the radiative decay rates can be particularly important for heavy mesons since, unlike the case of pions and kaons where such decays have been studied using Chiral Perturbation Theory (ChPT) [5][6][7][8][9], no model-independent calculations have ever been performed. Even in the case of light mesons, although the quoted ChPT calculations represent a first-principles approach to the problem, the low-energy constants entering in the final results at O(p 6 ) have been estimated in phenomenological analyses relying in part on model-dependent assumptions.
In Ref. [10] a strategy to compute QED radiative corrections to the P → ν (γ) decay rates at O(α em ) by starting from first-principles lattice calculations was proposed. The strategy has subsequently been applied in Refs. [11][12][13][14][15], within the RM123 approach [16,17], to provide the first non-perturbative model-independent calculation of the decay rates π − → µ −ν µ (γ) and K − → µ −ν µ (γ). In these calculations the real soft-photon contributions have been evaluated in the pointlike effective theory and, using the ChPT results quoted above, the SD corrections have been estimated to be negligible for these processes (see [10]). In the same phenomenological analysis it has been shown that the SD corrections might instead be relevant for the decays of pions and kaons into electrons. Moreover, by using the same single-pole dominance approximation as originally used in Ref. [18], SD contributions have been estimated to be phenomenologically important for decays of heavy-flavour mesons.
In this paper we present the first non-perturbative lattice calculation of the rates for the radiative decays P → νγ, where P is a pion, kaon, D or D s meson. We use the N f = 2 + 1 + 1 gauge ensembles generated by the European Twisted Mass Collaboration (ETMC) and analysed for mesonic observables in Ref. [19]. Preliminary results from this study were presented in Ref. [20]; the decays of bottom mesons will be studied in future papers. Note also that Kane et al. have presented preliminary results for the decays D + s → + νγ and K − → −ν γ, where ± represents the charged leptons and γ is a hard photon with energy in the range of about 0.5-1 GeV in Ref. [21].
The plan of the remainder of this paper is as follows. In Section II we introduce the basic quantities which enter in the amplitude for the leptonic decay of a pseudoscalar meson with the emission of a real photon; in particular we define the axial and vector form factors F A and F V .
We express the decay rates in terms of these quantities in Appendix A. In Section III we describe the general strategy that we followed to extract the amplitudes from suitable Euclidean correlation functions and discuss finite-time effects. The presence of discretisation effects which diverge at small photon momenta is demonstrated in Section IV and Appendix C, together with a strategy for subtracting them non-perturbatively. In Section V we present the numerical results for pions, kaons, D and D s mesons. Many formulae which are used in the paper are discussed and derived in Appendices A-C. Finally, in Appendix D we present some of our numerical results, including the correlation matrices, in a way which we hope may be useful to readers who wish to use them in phenomenological applications.

II. DEFINITION OF THE FORM FACTORS
The non-perturbative contribution to the radiative leptonic decay rate for the processes P → ν γ is encoded in the following hadronic matrix-element, see left panel of Fig. 1, H αr W (k, p) = r µ (k) H αµ W (k, p) = r µ (k) d 4 y e ik·y T 0|j α W (0)j µ em (y)|P (p) , where r µ (k) is the polarisation vector of the outgoing photon with four-momentum k, p is the momentum of the ingoing pseudoscalar meson of mass m P (p ≡ (E, p), E = m 2 P + p 2 and p 2 = m 2 P ). The operators are respectively the electromagnetic hadronic current and the hadronic weak current expressed in terms of the quark fields ψ f having electric charge q f in units of the charge of the positron; ψ U and ψ D indicate the fields of an up-type or a down-type quark and for the mesons considered in this study U can be either an up or a charm quark and D a down or a strange quark. In order to calculate the full amplitude one has also to consider the contribution in which the photon is emitted from the final-state charged lepton, see the right panel of Fig. 1. This latter contribution however, can be computed in perturbation theory using the meson's decay constant f P . Both contributions are included in the formulae for the decay rate given in appendix A.
The decomposition of H αr W (k, p) in terms of scalar form factors has been discussed in Ref. [10] (see also [22]). Here we adopt the same basis used in that paper to write The term in the last line of Eq. (3), which we write as H αµ pt (k, p), is the point-like infrared-divergent contribution. The other terms correspond to the so called SD contribution, H αµ SD (k, p). H αµ pt (k, p) saturates the Ward Identity (WI) satisfied by H αµ W (k, p) as explained in detail in Appendix C. The four form factors H 1,2 and F V,A are scalar functions of Lorentz invariants, m 2 P , p · k and k 2 . Eq. (3) is valid for generic (off-shell) values of the photon momentum and for generic choices of the polarisation vectors. The knowledge of the four form factors in the case of off-shell photons (k 2 = 0) gives access to the study of decays in which the pseudoscalar meson decays into four leptons. These processes are very interesting in the search of physics beyond the Standard Model and will be the subject of a future work. In this paper we concentrate on the case in which the photon is on-shell.
By setting k 2 = 0, at fixed meson mass, the form factors are functions of p · k only. Moreover, by choosing a physical basis for the polarisation vectors so that one has Once the decay constant f P and the two SD axial and vector form factors F A and F V are known, the radiative decay rate can be calculated by using the formulae given in appendix A. These formulae are expressed in terms of the convenient dimensionless variable where m is the mass of the outgoing lepton in the P → ν γ decay.
Our definition of the form factor F A differs from the definition, F B A , of refs. [21,23] We note that F B A includes the point-like infrared divergent contribution which totally dominates at low values of x γ thus obscuring the interesting structure-dependent contribution. For this reason we strongly advocate the use of our definition [10]. Moreover the sign of F V used in this paper is opposite to the one used in Ref. [21].

III. FORM FACTORS FROM EUCLIDEAN CORRELATION FUNCTIONS
In order to relate the hadronic matrix element to Euclidean correlation functions, the primary quantities computed in lattice calculations, it is useful to express the H αr W (k, p), defined in Eq. (1) in Minkowski space in terms of the contributions coming from the different time orderings. To this end, we define and perform the t y integral, whereĤ is the QCD Hamiltonian operator, E = m 2 P + p 2 is the energy of the decaying meson P and E γ = |k| is the energy of the outgoing real photon.
The important observation that makes the lattice calculation possible by using standard effective-mass/residue techniques is that the integrals over t y appearing in the definition of H αr W (k, p) can be Wick rotated to the Euclidean space without encountering any obstruction. Such obstructions arise whenever there are states propagating between the operators in the Tproducts that have energies smaller than the energy of the external states [24]. This doesn't happen in our case. For this reason H αr W,1,2 (k, p) can be rewritten in terms of Euclidean integrals, both of which are convergent for physical (non-vanishing) photon energies. In Eqs. (11) and below t y is a Euclidean time variable. Indeed, the hadronic state of lowest energy that can propagate between the two currents is the pseudoscalar meson with spatial momentum p − k (it appears in the time-ordering H αr W,1 ) and we have As a consequence, H αr W (k, p) can be rewritten as From this observation it follows that the hadronic matrix-element can be extracted from the Euclidean correlation functions where P = iψ D γ 5 ψ U is a Hermitian pseudoscalar interpolating operator having the flavour quantum numbers of the incoming meson. In Eq. (14), using the translational invariance of the correlation function, we have moved the origin in time to the pseudoscalar source, P (0, x), and placed the weak current at t.
In the large-t limit one has where the ellipsis represents the sub-leading exponentials.
The expressions for the correlation function C αr W (t; k, p) in Eq. (14) and for the ratio R αr W (t; k, p) in Eq. (15) refer to the ideal case of a lattice with infinite time-extent. The extraction of the matrix elements from correlation functions computed on a finite lattice in our numerical simulations is discussed in Appendix B. Although some of the details of the appendix refer to our specific lattice procedures (the choice of lattice Fermions, renormalisation of the operators, etc.) the strategy itself is general and can be directly translated to other lattice discretisations of QCD and of QED.

FIG. 2:
Schematic diagrams representing the correlation function C αr W (t, T /2; k, p) used to extract the form factors, see Appendix B. The interpolating operator for the meson P and the weak current j W are placed at fixed times 0 and t, and the electromagnetic current j em is inserted at t y which is integrated over 0 ≤ t y ≤ T , where T is the temporal extent of the lattice. The left and right panels correspond to the leading contributions to the correlation functions for t y < T /2 and t y > T /2 respectively, with mesons propagating with momenta p or p − k.
Here in the main text, we use Fig. 2 to illustrate the strategy used in our numerical simulations, performed with (anti-) periodic boundary conditions in time for the (fermionic) bosonic fields, to extract the form factors. The two panels in Fig. 2 represent the forward (0 t T /2) and backward (T /2 t T ) halves of the lattice. In both cases, the t y integral is dominated by the region in which t y is close to t, allowing for the propagation of the lightest state over the longest time interval.
In Fig. 3 we show two more diagrams to illustrate two important points concerning our numerical calculation of the correlation functions and of the form factors. The diagram in the left panel shows a quark-disconnected contribution to the correlation function originating from the possibility that the external real photon is emitted from sea quarks. In this work we have been using the so-called electroquenched approximation in which the sea-quarks are electrically neutral. In practice this means that we have neglected the contributions represented by the diagram in the left panel of we set the boundary conditions for the "spectator" quark such that ψ(x + nL) = exp(2πin · θ s /L)ψ(x), where L is the spatial extent of the lattice in each spatial direction. We treat the two propagators that are connected to the electromagnetic current as the results of the Wick contractions of two different fields having the same mass and electric charge but satisfying different boundary conditions [26]. This is possible at the price of accepting tiny violations of unitarity that are exponentially suppressed with the volume. By setting the boundary conditions as illustrated in the figure we have thus been able to choose arbitrary (non-quantised) values for the meson and photon spatial momenta by tuning the real three-vectors θ 0,t,s . We find that the most precise results are obtained with small values of |p| and in particular with p = 0.
The numerical results presented in the following sections have been obtained by setting the non-zero components of the spatial momenta along the third-direction, i.e.
With this particular choice of the kinematical configuration, a convenient basis for the polarisation vectors of the photon (see Appendix B for more details) is the one in which the two physical polarisation vectors are given by while the unphysical polarisation vectors vanish identically, 0 µ = 3 µ = 0. Notice that in this basis we have and, consequently, Using these formulae, we have built the following numerical estimators for the form factors, which we determine by fitting to the plateaux in the region 0 t T /2. The discussion here and below corresponds explicitly to the forward half of the lattice (0 t T /2).
We combine the results with those from the backward half (T /2 t T ) by exploiting timereversal symmetry as explained in Appendix B .
The ratios R jr W (t, T /2; k, p) appearing in Eqs. (21) and (22), which we evaluate separately for the axial (W = A) and vector (W = V ) components of the weak current, are the finite-T generalisations (see Eq. (B19)) of the ratios R αr W (t; k, p) defined above in Eq. (15). The values of the meson energies and of the matrix elements P |P |0 needed to build these estimators have been obtained from standard effective-mass/residue analyses of pseudoscalar-pseudoscalar two-point functions. We have also computed the pseudoscalar-axial two-point functions from which we have extracted the decay constants f P on our data sets in order to be able to separate the SD axial form factor F A from the point-like contribution 2f P /(m P x γ ).

DISCRETISATION EFFECTS
In this section we want to stress a very important issue associated with infrared divergent cutoff effects which can jeopardise the extraction of F A at small values of x γ . We also introduce a strategy to overcome this problem. Ref. [15].
In Fig. 4 we plot F A (x γ )+2f P /(m P x γ ), the sum of the point-like and SD axial form factors which is extracted directly from the correlation functions using R A (t) (see Eq. at small x γ . Using the decay constant and mass, computed in the standard way from the two-point functions, we can in principle subtract the point-like term and extract F A (x γ ). However, this turns out to be very difficult because of the possible presence of discretisation effects which cannot be excluded by the WI of the lattice action. Moreover, these lattice artefacts diverge as x γ → 0. We now propose a non-perturbative method to eliminate this problem.
At finite lattice spacing the axial form factor is constrained, as in the continuum (see Eq. (4)), by an exact lattice WI that is true at all orders in the lattice spacing a (see Appendix C). The label L here, and in the remainder of this section, stands for "Lattice" as the discussion concerns the Ward Identity in a discrete space-time. It should not be confused with the spatial extent of the Lattice. This however does not exclude the presence of cutoff effects in Eq. (21). These are terms of O(a 2 ) 1 and, in 1 We assume here that we are using a lattice discretisation in which the leading artefacts are O(a 2 ). For Wilson Fermions in which they are O(a), the discussion has to be modified accordingly.
where the ellipsis represents higher orders in a 2 , while the quantities ∆F A and ∆f P depend upon the parameters of the theory regularised on the lattice, on the light and heavy quark masses and upon Λ QCD . Discretisation effects in the pseudoscalar masses are also absorbed into ∆F A and ∆f P . The crucial point to notice is that the lattice decay constant of the WI in Eq. (23) f L P = f P + a 2 ∆f P . This implies the presence of the extra term of O(a 2 /x γ ) which appears, in spite of the naive expectations based on the exact lattice WI. Thus the coefficient of the last term in Eq. (24) is not in general given by 2f L P /(m P x γ ), where f L P is the quantity extracted from the axial-pseudoscalar lattice correlation functions at finite lattice spacing. More precisely, once the matrix element 0|j α A (0)|P (p) is parametrised as in Eq. (23), the definition of f L P at fixed cut-off depends upon the choice of the index α and of the lattice momentum p α L and, for this reason, is not unique. Therefore, given a generic definition of f L P , one cannot expect a complete cancellation of the infrared divergent term on the right-hand side of Eq. (24), because a residual lattice artefact will survive generating an effective, unphysical infrared divergent contribution to F sub A (x γ ) at finite cutoff (a 2 ∆f P = f P − f L P + a 2 ∆f P ). This phenomenon is illustrated for the D s meson in Fig. 5 where F sub A is plotted as function of x γ . Since the subtraction of the potentially divergent term is incomplete we observe a fast rise of the effective F sub A (x γ ) at small values of x γ . For this reason, even if one has data at different values of the lattice spacing, it is particularly difficult to extract the , especially at small x γ and for heavy mesons. This is illustrated by the intermediate (red) points in Fig. 5 which were obtained by fitting and subtracting the O(a 2 /x γ ) artefacts. The divergence at small x γ is reduced but the relative statistical uncertainties are increased.
We now present an alternative strategy that avoids this problem. In Appendix C we show that the correlation function C αr A (t; k, p) has a smooth behaviour as a function of k and that from C αr (20)). We can then construct 0 0.05 (25). The divergence at small x γ is reduced by fitting and subtracting the O(a 2 /x γ ) artefacts, at the price of increased uncertainties at small x γ ; these are the intermediate (red) points. The most accurate results are given by F N P sub A , obtained by the non-perturbative subtraction of these artefacts as in Eq. (27) and are shown by the lower (black) points. The data are obtained using the ensemble B55.32 of Ref. [15].

the quantityR
that, by construction, vanishes identically at x γ = 0. Up to statistical uncertainties, each term in the sums in the numerator and denominator of Eq. (26) is independent of the indices j, r. For the study of the constraints imposed by the electromagnetic Ward identity, it is helpful to view the right-hand side as H jr A (k, p)/H jr A (0, p) − 1 (which is also independent of j, r). From the improved estimatorR A (t) we can extract the structure dependent form factor F A using a quantity that we also show in Fig. 5 and that, in contrast to F sub A , does not show any divergent behaviour at small x γ . The reduction of the uncertainty on F A (x γ ) usingR A (t), with respect to a fit to the right-hand side of Eq. (25), as shown in Fig. 5, is impressive, particularly at small x γ and also for heavy mesons where there are discretisation effects of O(a 2 m 2 D (s) ). In the following we will only present results obtained with this method.
The knowledge of C jr A (t, T /2; 0, p) allows us also to define an alternative estimator for the form that we find has reduced statistical errors compared to R V (t). Note that because of parity symmetry the correlation function C jr V (t, T /2; 0, p) = 0, but this is only approximately true when it is estimated using a finite statistical sample. We find that taking the difference C jr in the numerator of Eq. (28) results in a significant reduction of the statistical uncertainty for physical values of x γ (see Eq. (7)).

V. NUMERICAL RESULTS
The results presented in this paper were obtained using the ETMC gauge ensembles with N f = 2+1+1 dynamical quarks at three different values of the lattice spacing, a = 0.0885(36), 0.0815 (30) and 0.0619(18) fm, with meson masses in the range 220-2110 MeV. Details about these ensembles are given in table II of Ref. [15], see also Table I   our simulation can also be found in Tables I and II in Appendix D (see also Table II of Ref. [15]).

Renormalisation of the corresponding axial-vector and vector currents with Twisted Mass Fermions
are the unrenormalised quantities as explained in Eq. (B6), Z A has been computed with method M2 and Z V with the WI [19]. In Table II  For pions and kaons, we have covered the full physical range of  (indeed we even have data for unphysical values corresponding to x γ > 1).
For the pion, guided by ChPT, we fit to the formula This is certainly not the most general formula to include higher orders terms in ChPT, for example it does not contain chiral logarithms, but it is sufficiently simple and adequate to describe the pion data. The two coefficientsc 0 andc 1 take into account possible mass-independent discretisation effects. c 1 is multiplied by m 2 π because it arises in higher orders in ChPT. On the other hand the discretisation term proportional toc 1 is not multiplied by the mass of the meson because at this order in a there is an explicit violation of chiral invariance in the lattice Fermion Lagrangian.
When using the simpler expression in Eq. (29) we exclude data at pion masses m π 350 MeV.
Since in our data we have pion masses up to about 500 MeV, we have also performed fits in the full range by modifying Eq. (29) to include higher order terms as follows: The higher-order coefficients ∆c 0 , ∆c 0 ,c 1 , ∆c 1 and ∆c 1 have very little effect on the extrapolated results and for this reason they are not well determined. Indeed they only contribute to a slight increase in the uncertainty in the value of the pion form factors at x γ = 0 and in the slope in x γ . Similarly, in the different fits that we performed, we also added some of the possible lattice artefacts that break Lorentz invariance, for example those proportional to a 2 |k| 2 (in the frame where the meson is at rest), where k is the momentum of the photon. We found that their effect is very small and this was only taken into account in the evaluation of the final uncertainties.
Since SU (3) breaking effects may be important, and we only have results obtained at two values of the strange quark mass, for the kaon we first interpolate the form factors to the physical kaon mass and then fit them to the formula with pion masses m π < 350 MeV and in the full range of pion masses. Formulae (31) and (32) for the kaon are equivalent to those in (29) and (30) respectively for the pion. The presence of the constant term c 1 in Eqs. (31) and (32) is a reflection of the fact that the strange quark mass is fixed to its physical value. To simplify the notation we have used the same symbols for the coefficients in Eqs. (29)-(32) but the reader should note that their values are different in each case. We do not have sufficient data to include terms proportional to m 2 K m 2 π or m 4 π with logarithmic corrections in Eq. (32). In Fig. 8 we present the values of the pion (left panels) and kaon (right panels) form factors F A (x γ ) (upper panels) and F V (x γ ) (lower panels) as a function of x γ for the configurations at a = 0.0619 fm. The plotted points with error bars correspond to different values of the light-quark mass at several values of x γ . The points with large uncertainties (σ F A,V ≥ 0.01 for the kaon or σ F A,V ≥ 0.008 for the pion) are shown with faint grey symbols. These points are obtained for mesons with substantial non-zero momenta p = 0. The results of our simulation are compared to the lowest order in ChPT, given by the results of the fits, using the formulae given in Eqs. (30) and (32), after extrapolation to the continuum limit and physical quark masses, together with the corresponding uncertainties.
where P represents π or K and we take (L r 9 + L r 10 ) 0.0017 [28]; this is indicated by the horizontal red lines. The blue lines and green bands are the results and uncertainties of the fits, obtained using Eqs. (30) and (32) after the extrapolation to physical quark masses and to zero lattice spacing has been performed. In Fig. 9 we show the value of the pion (left) and the kaon (right) form factors F A (x γ ) (upper) and F V (x γ ) (lower) as a function of x γ , extrapolated to the continuum at the physical point, either using Eqs. (29) and (31) for the pion and kaon respectively to fit the data, full green bands, or by using Eqs. (30) and (32) where the constants F 0 A,V are a function of the light quark masses. Since, however, for this exploratory study, we have only two values of the heavy quark mass, both around the charm mass, we prefer to interpolate the values of the form factors to the physical charm quark mass and then to fit the result with the simple formula We have also performed fits with the pole-like formula x γ at fixed lattice spacing (a = 0.0815 fm) for the ensemble B25.32 [15]. The full blue and shaded orange bands are the results of the fits with the polynomial or pole formulae given in Eqs. (35) and (36) respectively.
In this first study, we only have results for the D (s) mesons in the range 0 ≤ x γ ≤ 0.4, corresponding to E γ 400 MeV in the rest frame of the hadron. In Fig. 10 we give the results for the form factors of the D s meson, F A (x γ ) and F V (x γ ), at a = 0.0815 fm. The full blue and shaded orange bands are the results of the fits with the polynomial or pole formula given in Eqs. (35) and (36) respectively. Since the lattice spacing is fixed, the coefficientsd 0,1 are not included in the fit. We see that the both the fits give a good description of our results in the region where we have data, but differ significantly for x γ ≥ 0.4. This means that, although both the linear and the pole fits describe accurately the form factors in the region in which we have data, it is not reliable to use these fits in the region x γ ≥ 0.4. In our future investigations we plan to provide non-perturbative data for the form factors in the full kinematical range 0 ≤ x γ ≤ 1 − m 2 /m 2 D (s) . In Fig. 11  The polynomial and pole fits correspond to Eqs. (35) and (36) respectively.
with their central red lines are the results of a single fit to all the data after extrapolation to the continuum limit and to physical quark masses. The discretisation artefacts, which include ones of O(m 2 c a 2 ), while approximately of the expected size, appear to be relatively large because the form factors are small. In fact the form factors at the three lattice spacings we have at our disposal are fully consistent, within our uncertainties, with a linear behaviour in a 2 , as illustrated in Fig. 12 where the form factors at x γ = 0.2 are presented as a function of the lattice spacing. The points in the figure are obtained after extrapolation to physical quark masses either using a polynomial of pole ansatz corresponding to Eqs. (35) or (36) at fixed lattice spacing. In this first study, with only three lattice spacings at our disposal, we are unable to include corrections of higher order in a 2 beyond those present in Eqs. (35) and (36). In Appendix D we have estimated their effects in the uncertainties of our final results for the form factors.
We also study our physical results (i.e those obtained after the continuum and chiral extrapolations) as a function of x γ by fitting them to the following linear expressions: where P represents each of the pseudoscalar mesons, π, K, D and D s .
For the axial form factors we find: and for the vector form factors we obtain In Eqs. (38) and (39), for each of the C's and D's, ρ C,D is the correlation between them, defined by where C i and D i are the jackknife samples and the sum runs over all the jackknifes following the procedure in Appendix A of Ref. [29].
For the pion and kaon we can compare the constants C π,K A,V in Eqs. (38) and (39) with the constant (i.e. x γ -independent) values obtained in ChPT using Eq. (33): F π A = 0.0119, F π V = 0.0254, F K A = 0.042, F K V = 0.096.
In the remainder of this section we present a brief comparison of our results with experimental data. A more detailed phenomenological analysis will be presented in a separate paper.
For the kaon the PDG quotes the two combinations F K V ± F K A . They present separate values obtained from K → e decays, F K V +F K A = 0.133 (8), and from K → µ decays, F K V +F K A = 0.165 (13). Of course the results should be independent of whether the final-state charged lepton is an electron or muon. For this combination of form factors our value is F K V + F K A = 0.161 ± 0.013 at x γ = 0 and F K A + F K V = 0.1363 ± 0.0096 at x γ = 1. For the other combination of form factors the PDG quotes F K A − F K V = −0.21 (6) obtained from K → µ decays, which is quite different from our result The results in Eqs. (38) and (39) can be combined with the values of the decays constants computed in Ref. [19] and [30] f π = (130.41 ± 0.20) MeV to compute the differential or total decay rate using the expressions given in Appendix A.
For completeness, we also present the constantsC We found that the extraction of the axial form factor F A at small values of x γ is problematic because of the presence of very large discretisation effects of O a 2 / r 2 0 x γ and we provided a procedure for the non-perturbative cancellation of these systematic errors. We also found that for Simulations on one or more finer lattices would enable us to improve our estimates of the higher order artefacts and hence reduce the corresponding systematic uncertainty. Such preliminary studies of charmed mesons are also essential in order to study radiative decays of B mesons in the future. In this respect the use of the ratio method may also be very useful [31].
Although the present study clearly can and will be improved by, for example, increasing the statistics, covering the full range of x γ for D and D s mesons or simulating on a finer lattice, the results presented in this work already allow for an accurate comparison of the theoretical predictions with experimental measurements and we will discuss the phenomenological implications of our results in a forthcoming paper.
In future we also plan to study the emission of off-shell photons (k 2 = 0), computing all four form factors appearing in Eq. (3), which would allow us to predict the rates for processes in which the pseudoscalar meson decays into four leptons. These processes are very interesting in the search of physics beyond the Standard Model [32][33][34].
Appendix A: Expressions for the decay rates in terms of F V and F A In this appendix we present the explicit formulae needed to evaluate the total and differential decay rates at order α em , combining the non-perturbative determination of the virtual corrections computed with the approach of Ref. [10] with the calculation of the structure-dependent (SD) form factors F A and F V determined with the method proposed in this paper. These formulae can be used to compute the double differential decay rates d 2 Γ/(dx γ dx ), the single differential decay rates, dΓ/dx or dΓ/dx γ , as well as the integrated decay rate Γ(∆E γ ) = 2∆Eγ /m P 0 dx γ (dΓ/dx γ ) (∆E γ is the upper limit on the energy of the emitted photon in the meson rest-frame).
The exchange of a virtual photon depends on the hadron structure, since all momentum modes are included, and the amplitude must therefore be computed non-perturbatively. On the other hand, the non-perturbative evaluation of the amplitude for the emission of a real photon is not strictly necessary [10]. Indeed, it is possible to compute the amplitudes for real-photon emission in perturbation theory when x γ is sufficiently small that the internal structure of the decaying meson is not resolved. The infrared divergences in the non-perturbatively computed amplitude with the exchange of a virtual photon are cancelled in the decay rates by those present in the emission of a real photon, even when the latter is computed perturbatively. The reason for this cancellation is the universality of the infrared behaviour of the theory (i.e. the infrared divergences do not depend on the structure of the decaying hadron). For large photon energies, for example those present in the decays of heavy mesons, a full non-perturbative determination of the relevant amplitudes is necessary.
To calculate the partial rates for the emission of a hard real photon it is sufficient to know the SD form factors, F A and F V , and the meson's decay constant f P . For the integrated rate Γ(∆E γ ) instead, in the intermediate steps of the calculation it is necessary to introduce an infrared regulator.
To this end, in order to work with quantities that are finite when the infrared regulator is removed, it is very useful to organise the inclusive rate Γ(∆E γ ) = Γ(P − → −ν (γ))| Eγ ≤∆Eγ as follows 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 µ γ is an infrared regulator. On the right-hand side of Eq. (A1) the quantities Γ 0 (L) and Γ 1 (∆E γ ) are evaluated on the lattice.
The terms in the first parentheses on the right-hand side of Eq. (A1), Γ 0 (L) and Γ pt 0 (L), have the same infrared divergences which therefore cancel in the difference. Here we use the lattice size L as the intermediate infrared regulator by working in the QED L formulation of QED in a finite volume [35] but any other consistent formulation of QED on the lattice can also be used. The is independent of the regulator as this is removed [11]. Γ 0 (L) depends on the structure of the decaying meson and is computed non-perturbatively Refs. [11][12][13][14][15].
In the terms in the second parentheses on the right-hand side of Eq. (A1) the decaying meson is taken to be a point-like charged particle and both Γ pt 0 (µ γ ) and Γ pt 1 (∆E γ , µ γ ) can be computed directly in infinite volume, in perturbation theory, using some infrared regulator, for example a photon mass µ γ = m γ . Each term is infrared divergent, but the sum is convergent [2] and independent of the infrared regulator. In Refs. [10] and [11] the explicit perturbative calculations of Γ pt 0 (µ γ ) + Γ pt 1 (∆E γ , µ γ ) and Γ pt 0 (L) have been performed with a small photon mass µ γ or using the finite volume respectively, as the infrared regulators.
Finally, the term on second line of the right-hand side of Eq. (A1) is infrared finite. It can be computed in the infinite-volume limit requiring only knowledge of the structure dependent form factors, F A (x γ ) and F V (x γ ) and of the meson's decay constant f P where Γ SD is the structure-dependent contribution and Γ INT is that from the interference between the SD and point-like components of the amplitudes. Both Γ SD and Γ INT are separately infrared finite and there is no need to introduce an infrared regulator in this term.
We express the differential decay rate in terms of the following quantities: • the two dimensionless kinematical variables where m the mass of the lepton , with r = m /m P ; • the decay constant of the meson f P ; • the two SD axial and vector form factors F A and F V .
The differential decay rate is given by the sum of three contributions, where Γ (0) is the leptonic decay rate in the absence of electromagnetic corrections. This is given by where G F is the Fermi's constant and V CKM the relevant CKM matrix element.
The quantities in the braces on the right-hand side of Eq. (A4) are given by and correspond to the contribution of the point-like approximation, to the SD contribution and to the interference between point-like and SD terms respectively. The kinematical functions appearing in Eq. (A6) are given by The distribution with respect to the photon's momentum is obtained after integrating over the lepton's momentum As x γ → 0 the allowed kinematical range for x is squeezed around its maximum, x min (x γ ) = f pt (x γ , x ) ∼ 1/x 2 γ , all the other contributions vanish in the soft-photon region, which is consequently dominated by the pointlike (eikonal) result The 1/x γ behaviour of the differential rate at small x γ leads to a logarithmic infrared divergence in the total rate. It is cancelled by the infrared divergence in the O(α em ) virtual corrections to the inclusive decay rate. The SD and INT contributions vanish at small x γ .
Eqs. (A4)-(A9) allow us to compute the spectrum dΓ/dx γ . We advocate organising the determination of the integrated rate in terms of the three sets of parentheses on the right-hand side of Eq. (A1). The procedure to evaluate the term in the first parentheses, Γ 0 (L) − Γ pt 0 (L), is explained in detail in Ref. [10], where the explicit expression for the term in the second parentheses, Γ pt 0 (µ γ ) + Γ pt 1 (∆E γ , µ γ ), can also be found. The third term on the right-hand side of Eq. (A1), is the subject of this paper. As explained above, both Γ SD and Γ INT are infrared finite and are obtained by integrating the differential rates over the physical range of x γ ,

Appendix B: Calculating matrix elements from finite Euclidean lattices
In this appendix we derive some useful formulae for the extraction of the two relevant form factors, F A,V , from the Euclidean correlation functions expressed in terms of lattice operators on a lattice with finite time extent T .
In order to construct the finite T equivalent of C αr W (t; k, p) in Eq. (14), that we will denote as C αr W (t, T /2; k, p), it is convenient to define the following hadronic correlation function at fixed t and where . . . LT denotes the average over the gauge field configurations at finite L and T and we introduced suitable independent vectors r (k), r = 1, 2, corresponding to the physical polarisations of the emitted photon. A possible simple choice, and one in which the unphysical polarisations vanish explicitly, is given by The polarisation vectors satisfy Since in our simulations we always use k = (0, 0, |k|), the polarisation vectors reduce to The "topology" of the correlation function in Eq. (B1) is explained in Figure 2: • the incoming meson is interpolated at fixed spatial momentum p by the pseudoscalar operator P placed at time t = 0 • the hadronic weak current j α W (t) is placed at the generic time t. We used a local discretisation of the weak current that, in the Twisted-Mass discretisation of the fermionic action used in this work [36], is explicitly given by where j α V (t) and j α A (t) are the vector and axial components that include the corresponding renormalisation factors. Note that the renormalisation factors to be used in Twisted-Mass at maximal twist are chirally-rotated with respect to the ones of standard Wilson fermions [37].
In Eq. (B6) ψ U indicates the field of an up-type quark that, for the mesons considered in this study, can be either an up or a charm. Similarly, ψ D can be either a down or a strange quark field. The actions of the up-type and down-type quark fields have been discretised with opposite values of the chirally-rotated Wilson term in order to numerically suppress O(a 2 ) lattice artefacts in the meson masses [37,38]; • the electromagnetic current j µ em (t y , y), carrying a three-momentum k is inserted at y = (t y , y). This current is defined by where f is the flavour index, q f is equal to 2/3 for up-type quarks and to −1/3 for down-type quarks. and In Eq. (B8) U µ (x) are the QCD link variables and the signs ± are induced by the choice made in the case of the flavour f for the sign of the chirally-rotated Wilson term [17].
We have used e −ik·(y+î/2) , rather than the simpler, standard exponent e −ik·y , for the Fourier transform of the current appearing in Eq. (B1) Our choice of the exponent, which is equivalent to standard one in the continuum limit, is more convenient for the discussion of the lattice WIs since we have used the point-split exactly conserved electromagnetic current in our simulations.
• A technical subtlety needs to be stressed here. As discussed in the main text, in order to choose arbitrary (non-discretised) values of the spatial momenta for the meson and for the photon, we have introduced a "flavoured" extension of the electromagnetic current (see the explanation in the caption of Figure 3). In practice, in order to have two quarks (ψ 0 and ψ t , where 0 and t are labels for the quark fields) having the same mass, the same electric charge, the same sign of the chirally-rotated Wilson term but different boundary conditions [26], the expression to be used in the numerical calculation is where we have used the fact that (see Eq. (17)) and we have defined, as usually done in implementing twisted boundary conditions [25], the periodic fields In all the formulae that will follow the range of the time parameters is extended over the full lattice extension, 0 ≤ t < T and 0 ≤ t y < T .
We are now ready to define the finite T correlation function In the continuum and large-T limits one can readily show that for 0 t T /2 where C αr W (t; k, p) is the correlation function introduced in the main text defined in Eq. (14), H αr W (k, p) is the physical matrix element defined in Eq. (1) and the ellipsis represent sub-leading exponentials.
For negative time t on the other hand, i.e. for time separations such that T /2 t T , in the continuum and large-T limits we have with the ellipsis again representing the sub-leading exponentials.
It is useful to note that, in order to separate the axial and vector form factors, it is enough to compute separately the correlation functions corresponding to the vector, C αr V (t, T /2; k, p), and the axial, C αr A (t, T /2; k, p), components of the weak current. Moreover, from the properties we deduce the following properties of the corresponding correlation functions under time reversal: Using these relations, the quantities were extracted from the ratios of the correlation functions averaged over the two temporal halves of the lattice In all the formulae of this Appendix we have used continuum notation for the four vectors but the momentum and energy carried by the current (including the associated projectors) have to be read by performing the following substitutions In the lattice regularisation that we are using (i.e. Wilson quarks at maximal twist), Eqs. (B16) and (B17) hold for given values of the indices α and r ∈ {1, 2} only up to O(a 2n+1 ) lattice artefacts (for integer n). One can show however, that, as a consequence of exact lattice symmetries (see e.g. Refs. [37,38]) and the choice of momenta and polarisation vectors given in Eqs. (17) and (18) which are precisely those that occur in Eqs. (26) and (28)   In this appendix we study the Ward Identity (WI) that relates the axial correlation function C αr A (t; k, p) to the axial-pseudoscalar correlation function and, consequently, the matrix element H αr A (k, p) to the decay constant of the meson f P . As discussed in the main text, a careful analysis of the cut-off effects reveals that the WI does not exclude the possibility of different O(a 2 ) artefacts appearing in the decay constant extracted from the three-point function and that from the twopoint function.
We start with a remark about the matrix element of the axial current, determined at a finite lattice spacing a and using a particular lattice discretisation, which we write in the form where the combination f L P p α L is a vector under the orthogonal group H(4) 2 . At finite a the definition of the lattice decay constant f L P depends on the definition that we assume for the lattice momentum p α L , for example we may choose p α L = p α or p α L = 2/a sin[ap α /2], where p α is the continuum value of the momentum. In particular we definef P by 0|j α A (0)|P (p) =f P (p)p α . Note thatf P = f P + O(a 2 ), where f P is the continuum value of the decay constant.
Consider the following correlation function which is relevant to our study, With Wilson-like Fermions, such as those used in our study, at fixed lattice spacing (for simplicity in the T → ∞ limit), the electromagnetic WI implies that [39] 1 a where integrals over the spatial coordinates have to be read as lattice sums and, in the case of a real photon, k 0 = iE γ (k) = i|k|.
The WI can be rewritten in the form where we have defined (note the shift in the exponent with respect to Eq. (C3)) and We can derive the Ward identity for the matrix element itself by going onto the mass shell of the pseudoscalar meson, which in the Euclidean corresponds to selecting the energy of the external hadronic state to be E P as | − t| becomes very large.
Consider first the case with k = 0. In this case the second term on the right-hand side of Eq. (C4) does not contribute because it corresponds to a different energy. Thus, in this case, we have the following identity which is true at all orders in a: Note that to arrive at this identity we do not need to specify the choice of f L P . As a → 0 the discretised derivative in Fourier space 2/a sin(ak µ /2) → k µ and we recover the continuum WI in where ∆ −1 is some version of a lattice boson propagator, for example as a → 0, A(k, p) = 1 + O(a 2 ) and T αµ (k, p) = (2p − k) µ (p − k) α + O(a 2 ) are functions of the momenta which depend on the lattice regularisation and f L P is the meson decay constant extracted from the matrix element in Eq. (C1). We therefore have At fixed lattice spacing the only condition that must be satisfied is that in applying the WI, the denominator ∆ of Eq. (C9) disappears. The WI guarantees that the right-hand side of Eq. (C12) is the matrix element 0|j α A (0) |P (p) including all orders in a. By iterating order by order in a, we may find solutions of the form From the above discussion we conclude that the lattice H jr A (k, p) has the form where the dots represent higher order discretisation corrections.
In order to implement the strategy described in Eq. (27), we need to perform a direct calculation of H jr A (0, p) and hence to study the k → 0 limit of the WI. The problem is non-trivial because from the spectral analysis of C αµ A (t; k, p) it follows that where the ellipsis represent exponentially suppressed contributions, with a gap that is O(m π ). The first exponential corresponds to the on-shell external meson P and represents the state we are interested in. The second exponential corresponds to the state P + γ where both the meson and the photon are on-shell and have a total momentum p and a relative momentum k. A similar time dependence also appears in the WI from the rotation of the P source; this is the second term on the right-hand side of Eq. (C4). As discussed above, when k = 0 it is possible to isolate the matrix element corresponding to the state P .
The problem we now address is to study the limit k → 0, paying special attention to the leading cut-off effects. This can be done by using the exact WI satisfied by C αµ A (t; k, p) at finite lattice spacing; in particular we aim to understand the structure of the correlation function C αµ A (t; k, p) at k = 0. To this end we consider the two-point correlation functions on the right-hand side of Eq. (C4) when α is a spatial index (the case α = 0 is similar, but in the following we shall concentrate on the case α = 1, 2, 3). By setting E γ = 0 in the last term of Eq. (C6), we have where the ellipsis represent sub-leading exponentials (the gap is at least 2m π ). In the previous expressionsf P (p) = f P + O(a 2 ),Ĝ P (p) = G P + O(a 2 ) andÊ P (p) = E P (p) + O(a 2 ) where f P , G P and E P (p) are respectively the continuum decay constant, the continuum matrix element of the pseudoscalar density used as interpolating operator, G P = 0|P |P (p) , and the continuum energy of the meson.
By using the previous two expressions and by differentiating Eq. (C4) with respect to the component k i of k and then setting k = 0 we obtain where the ellipsis again represents the sub-leading exponentials and we have used the fact that As can be seen, the structure of C ji A (t; 0, p) is highly non trivial. Note in particular the term linear in t that is a manifestation of the singular behaviour at large distances of the correlation function (this generates a double pole in momentum space that is at the original the infrared divergence). An important consequence of the strategy proposed in section IV which we have used in our calculations, is that the terms in squared brackets in Eq. (C17) disappear at any value of p when we contract the correlation function with the physical polarisation vectors of the photon.
With our choice of kinematics, these satisfy the relation Indeed, the H(3) symmetry implies that and thus We conclude that C jr A (t; 0, p) can be analyzed as expected to extract the coefficient of the leading exponential. We stress that the above demonstration shows that from C jr A (t; 0, p) we can extract precisely the decay constant which appears in the lattice matrix element of the axial current in Eq. (C1), without to have to make a choice for the lattice momentum p α L .

Appendix D
In this Appendix we present some numerical information that may be useful to the reader. We start by listing in Tables I and II the Table II. the results obtained using Eqs. (37) and (42)