Pion and Kaon Distribution Amplitudes from Lattice QCD

We present the state-of-the-art lattice QCD calculation of the pion and kaon light-cone distribution amplitudes (DAs) using large-momentum effective theory. The calculation is done at three lattice spacings $a\approx\{0.06,0.09,0.12\}$ fm and physical pion and kaon masses, with the meson momenta $P_z = \{1.29,1.72,2.15\}$ GeV. The result is non-perturbatively renormalized in a recently proposed hybrid scheme with self renormalization, and extrapolated to the continuum as well as the infinite momentum limit. We find a significant deviation of the pion and kaon DAs from the asymptotic form, and a large $SU(3)$ flavor breaking effect in the kaon DA.

Introduction: Light pseudoscalar mesons play a fundamental role in quantum chromodynamics (QCD) as they are the (pseudo) Nambu-Goldstone bosons associated with dynamical chiral symmetry breaking (DCSB) [1], an important non-perturbative phenomena in the standard model (SM). Their internal structure and its impact on experimental measurements have been actively investigated for many years.
Among others, the leading-twist pion and kaon lightcone distribution amplitudes (DAs) are the simplest physical quantities to describe such internal structure, and provide a probability amplitude interpretation on how the longitudinal momentum of the pion/kaon is distributed among quarks in its leading Fock state. They are critical inputs for the description of hard exclusive reactions, such as the B meson weak decays [2,3] which provide useful information on CP violation and the Cabibbo-Kobayashi-Maskawa matrix, and play a crucial role for probes of new physics [4]; they are also important for the study of the pion elastic form factors [5], the pion-photon transition form factor [6], and of hard exclusive meson production which may give access to nucleon generalized parton distributions [7,8].
In the high energy limit, it is well-known that these DAs follow a simple asymptotic form φ(x) = 6x(1−x) [9]. However, their shapes at lower scales have been long debated. For example, the CZ model [10] proposed a "double-humped" shape for the pion DA which allowed a consistent description of the experimental data at that time, while other models (see, e.g., [11][12][13][14][15]) do not favor such a structure. In addition, the measured electromagnetic form factors of pion/kaon at relatively large momentum transfer also exhibit some puzzling feature that might be connected to their non-asymptotic behavior [16,17]. Theoretical calculations from lattice QCD will be able to shed more light on the shape of these DAs.
There have been a lot of lattice studies on the pion/kaon DA using the traditional moments approach [18][19][20][21][22][23][24]. The recently proposed large-momentum effective theory (LaMET) [25][26][27] allows to access the entire x-dependence of the DAs from first principle lattice calculations, instead of only the first few moments (for other proposals with applications to the DAs, see Refs. [28][29][30][31]). Using LaMET, several calculations of the x-dependence of meson DAs have been carried out [32][33][34]. However, a recent analysis [35] showed that the nonperturbative renormalization of the quasi-light-front (quasi-LF) correlation in LaMET could be highly nontrivial, especially when off-shell quark matrix elements are used. In such a case, even after renormalization there may still be residual linear divergences rendering the continuum extrapolation problematic. To resolve this issue, a self renormalization strategy [36] has been proposed, where one fits the divergence structure to a quasi-LF correlation and uses it for renormalization. The present work provides the first full implementation of this strategy, and shows that it indeed gives promising results.
Lattice simulation: The leading-twist light-cone DA of a pseudoscalar meson is defined as where U (0, ξ − ) = P exp ig s 0 ξ− ds n · A(sn) is the pathordered gauge link defined along the minus light-cone direction n (n 2 = 0). To extract this quantity, we calculate the following quasi-LF correlation on the lattice with momentum P along the z-direction [32]: where O Γ1 (z; y, t) ≡ψ 1 ( y, t)Γ 1 U ( y, y−zẑ)ψ 2 ( y−zẑ, t) is the quasi-LF operator withẑ being the unit vector in the z-direction, U ( x, x− z) is the spatial Wilson line connecting lattice sites x and x − z, ψ 2 Γ 2 ψ 1 is the interpolating field of the meson m, and Γ 1,2 are chosen as Γ 1 = γ z γ 5 , Γ 2 = γ 5 for the pseudoscalar meson. The ground-state matrix elements can be extracted from the following two-state fit formula: where H B m (z) is the normalized ground-state matrix element, c m and ∆E are free parameters accounting for (one or more) excited state contamination, which are exponentially suppressed in the large time limit. Based on the comparison of one-and two-state fits (see supplemental material [37]), we use the one-state fit results in the analysis below with t min = 0.72, 0.54, 0.42 fm (for P z = 1.29, 1.72, 2.15 GeV) which is large enough to eliminate the excited states contamination.
In this work, the simulation is done using the clover fermion action on three ensembles with 2+1+1 flavors of highly improved staggered quarks (HISQ) generated by the MILC collaboration [38,39], at physical pion mass with three lattice spacings 0.057, 0.088 and 0.121 fm. Hypercubic (HYP) smeared fat links [40] are used in both the fermion action and the quasi-LF operators in C m 2 to improve the signal-to-noise ratio. The rest of the simulation setup is summarized in Table I. In addition, we use momentum smeared 2-2-2 grid sources, repeat the calculation at several time slices, and average the forward and backward correlation functions to improve statistics. In total, we have 570 (cfg.) × 8 (grid source) × 8 (source time slices) × 2 (forward/backward), 730 × 8 × 6×2 and 970×8×4×2 measurements at three ensembles with a = 0.057, 0.088 and 0.121 fm, respectively.
Hybrid scheme with self renormalization: The bare quasi-LF correlation calculated above contains both linear and logarithmic ultraviolet (UV) divergences which have to be removed by renormalization. On the lattice, the numerical subtraction of linear divergences is extremely delicate. In particular, such divergences may not be fully removed if the RI/MOM renormalization scheme is used [35]. Here we adopt the self renormalization proposed in Ref. [36], which amounts to fitting the bare quasi-LF correlation and subtracting the relevant UV divergences. To be more precise, one fits the bare quasi-LF correlation at given hadron momentum and multiple lattice spacings with a perturbative-QCD-dictated parametrization that contains a linear divergence, a logarithmic divergence, and discretization effects. After removing all the UV divergences and discretization effects, one is left with the renormalized quasi-LF correlation encoding the intrinsic non-perturbative physics.
As suggested in Ref. [36], the UV divergences in the quasi-LF correlator can be determined by using, e.g., the pion PDF matrix elements M(z) ≡ π|O γt |π in the rest frame at multiple lattice spacings, and fitting the bare data M B to the following form [36] with the renormalization factor parametrized as [36] where the first term in the curly bracket is the linear divergence, m 0 denotes a finite mass contribution arising from renormalon ambiguity, etc., and f (z)a accounts for the discretization effects (The O(a) correction here arises from the mixed action effect in using clover valence fermions on HISQ sea ones). The last two terms come from the resummation of leading and subleading logarithmic divergences, which only affect the overall normalization at different lattice spacings. To partially account for higher-order perturbative effects as well as remaining lattice artifacts, we also treat d and Λ QCD as fitting parameters [36]. The renormalized matrix element is required to be equal to the continuum perturbative MS result at short distances (chosen to be z ∈ z s = [0.06, 0.18] fm as defined in [41]), which helps the determination of m 0 and d. In the calculation we use the MS renormalization scale µ = 2 GeV and Λ MS = 0.24 GeV.
In the present case, we follow the same strategy as above, except that the renormalized matrix element in the MS scheme is now required to be matched to the continuum perturbative MS result of the normalized quasi-DA matrix element at short distances in the rest frame, which reads at one-loop Z self turns out to be the same as Z self except for the value of d.
In Fig. 1, we show a comparison between the self renormalized quasi-LF correlations Re[e It is worth pointing out that the self renormalization strategy above does not apply at very small z due to finite lattice spacing artifacts in the data. In the ratio scheme [42], some degree of cancellation happens in the bare correlations between large momentum states and non-perturbative lattice renormalization factors. However, in the present case, the agreement of the selfrenormalized LF-correlation with the perturbative result extends down to z ∼ 0.06 fm which is our smallest lattice spacing. Thus, we only need to supplement it with the renormalized quasi-LF correlation at z = 0 which is normalized to 1. In this way, we obtain the fully renormalized quasi-LF correlation. To facilitate the subsequent matching procedure, we define a modified renormalized correlation by further dividing out the perturbative factor H MS, 1−loop m (z) so that the ratio scheme matching applies, Note that this is equivalent to using the hybrid renormalized quasi-LF correlation and matching, as the perturbative difference in the quasi-LF correlation is exactly compensated by that in the matching. From Fig. 1, we can see that the uncertainty of the renormalized quasi-LF correlation grows rapidly at large distance. A brute-force truncation of the correlation introduces unphysical oscillations [34] in momentum space after Fourier transformation. To resolve this issue, we adopt a physics-based extrapolation form [41] at large quasi-LF distance (λ = zP z ): where the algebraic terms in the square bracket account for a power law behavior of the DAs in the endpoint region and the exponential term comes from the expectation that at finite momentum ( P ) the correlation function has a finite correlation length (denoted as λ 0 ), which becomes infinite when the momentum goes to infinity. In this work, the Lorentz boost factor γ for the pion at the physical point is very large {9.21, 12.29, 15.36}, and thus the correlation length is very large. We therefore drop the e −λ/λ0 factor, and directly perform a polynomial extrapolation as suggested in [41]. The details of this extrapolation can be found in the supplemental material [37]. Numerical results: In Fig. 2 we show, as an example, the imaginary part of e izPz/2 H R m (z) for the pion (upper panel) and kaon (lower panel) at different lattice spacings with P z = 2.15 GeV. It reflects the SU(3) flavor breaking effects between the valence quarks in the light meson. For the pion it is consistent with zero within errors as expected, since we have used degenerate valence u/d quark masses in the ensembles. While in the case of kaon there is a non-vanishing imaginary part. Such an imaginary part increases slightly with P z , as observed in previous DA studies using LaMET [34,41], and a comparison of the results at different momenta can be found in the supplemental material [37].
The factorization can be done either in momentum space [43,44] or in coordinate space. Here we choose the latter, which results in where we take renormalization scale and factorization scale to be the same and set µ = µ R = 2 GeV in this paper. h R m is the LF correlation related to the light-cone DA through the following Fourier transformation h R m (x, y, λ, µ) = The perturbative matching kernel C up to the next-toleading order is given in the supplemental material [37].
The impact of the perturbative matching is illustrated in Fig. 3, where a Fourier transformation to momentum space has been performed. As can be seen from the figure, the matching broadens the quasi-DA in the physical region. Outside the physical region (x < 0 or x > 1), there still exists a non-vanishing tail, indicating potential effects of higher-order matching and higher-twist contributions. Nevertheless, in the unphysical region, the results are consistent with zero within ∼ 2 standard deviations.
With the results for P z = 1.29, 1.72, 2.15 GeV above, we can perform an extrapolation to P z → ∞ using the functional form: The final results of the π, K DAs are given in Fig. 4, where systematic uncertainties from renormalization scale, algebraic extrapolation, continuum and infinite momentum extrapolation have been taken into account. As the endpoint regions can not be reliably predicted by LaMET, we adopt a phenomenological x a (1−x) b extrapolation in this region (taken as 0 < x < 0.1 & 0.9 < x < 1). For comparison, we also plot the asymptotic form 6x(1 − x) and results from QCD sum rules [45], Dyson-Schwinger equations (DSE) [46], and reconstructed from moments calculations (OPE) [24]. As can be seen from the figure, both π and K DAs deviate significantly from the asymptotic form, but are close to the results from DSE and OPE calculations. The shape of π DA is much broader than the asymptotic form which could be explained as a direct expression of DCSB, as discussed in [14]. The relatively large uncertainties in Fig. 4 are expected to be reduced if the analysis of higher-order and higher-twist effects becomes available in the future. Summary: We present a state-of-the-art lattice calculation of π and K DAs using LaMET. The renormalization is done in the hybrid scheme with self renormalization proposed recently. Based on the results at physical light and strange quark masses with three lattice spacings and momenta, we perform an extrapolation to the continuum and infinite momentum limit. The final re-sults exhibit a significant deviation from the asymptotic form, while they are close to the DSE and OPE results, especially in the middle x region where our method is reliable. However, there are still some significant differences in the endpoint regions. This could be due to missing higher-power/ high-order corrections in LaMET which can be improved in future calculations, or due to effects of higher moments ignored in the OPE and DSE calculations.

ACKNOWLEDGEMENT
The calculations were performed using the Chroma software suite [47] with QUDA [48][49][50]  Additional information on the data analyses Dispersion relation: The effective mass can be extracted by the two-point correlation function at z=0 through a two-state fit. The dispersion relation of π and K can be derived from the parametrized form : Here C 3,π = −0.178 (19), C 3,K = −0.151 (13) parameterize the discretization effects. C 2,π = 0.9921(61), C 2,K = 1.0210(62) are consistent with the speed of light within 3σ. Excited state contamination: The common approach to extract the ground-state matrix elements H B m (z) is the two-state fit in Eq.(2) of the main text. However, to correct for excited state contaminations in this way requires high precision. An alternative solution is to use a one-state fit. As the excited state contaminations get suppressed when t becomes large, we can eliminate them by using a large enough t range.
As shown in Fig. 5, the excited state introduces a slight bend to the two-state fitted curve, especially in the small t region. However, with current accuracy, the excited state contaminations cannot be fitted effectively, and the one-state fit result is consistent with the two-state fit ones within statistical uncertainties. Large λ extrapolation: In the hybrid scheme, we adopt a physics-based extrapolation [41] at large quasi-LF distance (λ = zP z ). Depending on the hadron momentum, two extrapolation forms have been proposed: where the exponential term comes from the expectation that at finite momentum the correlation function has a finite correlation length (denoted as λ 0 ), which goes to infinity when the momentum goes to infinity. In this work, for light π, K meson at the physical point, the Lorentz boost factor is very large, so is the correlation length. We therefore use the formula in the second row in Eq. (14) for the extrapolation. For the pion, one can choose c 1 = c 2 , a = b due to isospin symmetry, and fit them using available lattice results in the large λ region. The pion result at P z = 2.15 GeV is shown in Fig. 7. We can see that the polynomial fit (blue/green diamonds) agrees with the original lattice data (red dots) well. The situation is similar in the case of the kaon. To estimate the modification effects of extrapolation form, we take two extrapolation regions(λ L1 > 8.2, λ L2 > 11 in Fig. 7 case) and treated their difference as an estimate of systematic error from extrapolation. Renormaliation scheme dependence: In this work, we have used a hybrid renormalization scheme based on self renormalization, where the renormalization factors are determined by fitting to the bare hadron matrix element. In a previous work [52], a hybrid scheme based on RI/MOM was used, where the fitting was done for the off-shell quark matrix element. As the lattice spacings used there are not very small, the linear divergences appear to be canceled up to numerical uncertainty. However, the fitting form is not well justified theoretically. Nevertheless, here we make a comparison between the two schemes, and show the results in Fig. 8. As can be seen from the figure, the renormalized matrix elements in the two schemes deviate less than 1σ in most regions. Thus, the scheme from [52] still yields consistent results with the hybrid scheme based on self renormalization, although the latter is theoretically more self-consistent. For comparison purposes, we also show the original renormalized π correlation without phase rotation in Fig. 10.
P z dependence of Kaon DA in momentum space: The momentum dependence of the continuum extrapolated results is shown in Fig. 11 for the kaon. Since the endpoint region behavior cannot be reliably predicted by LaMET, we shade the regions x < 0.1 and x > 0.9. As shown in the plot, the kaon DA is asymmetric around x = 1/2 with the strange quark having on average a larger momentum fraction (x is the momentum carried by the u/d quark). This asymmetry increases slightly with P z in agreement with our expectation, since s quark is heavier than u/d quark. A similar asymmetry has also been observed in previous DA studies in LaMET [34,41].
Prediction of the moments: In phenomenological research, the DAs are usually expanded in terms of the Gegenbauer polynomials: where a n s are Gegenbauer moments and C n s are Gegenbauer polynomials. One can extract Gegenbauer moments from our results with this equation: a n = 2(2n + 3) 3(n + 1)(n + 2) (16) and the moments up to a 4 are given in Table. II. The odd moments of π are identically zero due to the chargeconjugation invariance. The second moments from our calculation are ξ π 2 = 0.300(41), ξ K 2 = 0.258(32) (ξ n = 1 0 dx(2x − 1) n φ(x)), which agree better with a previous OPE calculation [19] rather than the recent one [24].

Matching in Coordinate Space
The matching in coordinate space in the MS scheme can be expressed as: with C x, y, where µ is the factorization scale, x, y are the Feynman integral parameters, z is the non-local separation, λ = zP z is the quasi-LF distance for the quasi-DA. h m is the LF correlation related to the light-cone DA φ(u, µ) through the following Fourier transformation (h m here has a different convention from h R m in Eq. (10) of main text) Such a matching is related to the matching in momentum space [43,44] by Fourier transformation.
In practice, we extract the DA by applying an inverse matching on the quasi-DA in the ratio scheme, with the matching kernel being given by C ratio x, y, z 2 , µ = C x, y, z 2 , µ After some manipulations, the final matching in coordinate space can be expressed as: where λ is the LF distance for the DA and:  12 shows a comparison between the coordinate and momentum space matching implementations. The results agree within errors in the physical region. The slight difference is due to different extrapolation choices, as in the former case we can do the large λ extrapolation after the matching, while in the latter the large λ extrapolation has to be done before Fourier transform and matching. The results from momentum space matching suffer from an oscillation at the endpoints which causes a more significant tail in the unphysical region. Such a behavior comes from artifacts related to the discontinuity of the momentum space matching kernel at the endpoints.