Light-cone distribution amplitudes of pseudoscalar mesons from lattice QCD

We present the first lattice determination of the two lowest Gegenbauer moments of the leading-twist pion and kaon light-cone distribution amplitudes with full control of all errors: $a_2^\pi=0.101^{+24}_{-24}$ for the pion; $a_1^K=0.0533^{+34}_{-35}$ and $a_2^K=0.090^{+19}_{-20}$ for the kaon. The calculation is carried out on 35 different CLS ensembles with $N_f=2+1$ flavors of dynamical Wilson-clover fermions. These cover a multitude of pion and kaon mass combinations (including the physical point) and 5 different lattice spacings down to $a=0.039\,$fm. The momentum smearing technique and a new operator basis are employed to reduce statistical fluctuations and to improve the overlap with the ground states. The results are obtained from a combined chiral and continuum limit extrapolation that includes three separate trajectories in the quark mass plane.


Introduction
Hadron light-cone distribution amplitudes (LCDAs) have been introduced four decades ago [1][2][3][4][5][6][7] in the context of the QCD description of hard exclusive reactions. The LCDAs are scale-dependent nonperturbative functions that can be interpreted as quantum-mechanical amplitudes. Within this article we will use the term "LCDAs" synonymous with the leading-twist LCDAs. The latter describe the distribution of the longitudinal momentum amongst the quarks in the leading Fock state contribution of a hadron wave function at small transverse parton separations. The pion LCDA is both the simplest LCDA and also the most important one in phenomenological applications. Unsurprisingly, it has received the most attention in the literature. Its precise knowledge is becoming increasingly relevant in flavor physics (where weak decays, such as B → π ν , B → ππ, etc., are providing information on the Cabibbo-Kobayashi-Maskawa matrix), in two-photon hard reactions (like γ * → γπ or γγ → ππ), and -as a tool to access the flavor separation in the nucleon generalized parton distributions -in hard exclusive electro-production (eN → eN π) with Bjorken kinematics. φ π (x) Figure 1. Models for the pion LCDA: the blue line shows the asymptotic shape corresponding to the limit µ → ∞, while the orange line depicts the CZ model [8,10] for µ = 0.5 GeV. The green point shows the QCD light-cone sum rule result [16] for the mid-point at µ = 1 GeV.
Theoretical attempts to predict the shape of the pion LCDA φ π (x, µ 2 ) as a function of the longitudinal momentum fraction x at a scale µ have a long history. The discussion was shaped for many years by the famous paper by Chernyak and Zhitnitsky (CZ) [8] who calculated the second moment in x of the pion LCDA using QCD sum rules [9] and found a number much larger than the result expected at asymptotically large scales. Based on this calculation, CZ proposed a particular model for the pion LCDA at low scales, known as the CZ model. Assuming the validity of perturbative QCD factorization, this model allowed for a consistent description of all experimental data on hard exclusive processes that were available at that time [10]. In figure 1 we compare the asymptotic LCDA φ π (x, µ 2 ) µ→∞ −−−→ 6x(1 − x) [4,5] with the CZ model. The latter corresponds to a doublepeaked distribution, where one of the constituents is most likely to carry a small (∼ 0.15) and the other one a large (∼ 0.85) fraction of the longitudinal pion momentum.
The CZ model received some criticism. On the one hand, the validity of collinear factorization in hard exclusive reactions at relatively low momentum transfer was questioned [11,12] and the role of a competing "soft" or "end-point" mechanism was emphasized. In particular it was shown [12,13] that the data on the pion form factor at Q 2 ∼ 1-3 GeV 2 could be described by the soft contribution alone, without any "hard" corrections. On the other hand, it was argued that the QCD sum rules employed in ref. [8] were not reliable as they may suffer from large contributions from operators of higher dimension. A model for such higher-order contributions using the concept of nonlocal vacuum condensates [14] yielded a much smaller value of the second moment than the CZ model, see ref. [15] for a state-of-the-art study. Finally, the explicit calculation [16] of the value of the pion LCDA at the mid-point x = 1 2 , using an at that time novel method, the light-cone sum rule (LCSR) technique, gave a rather large number, see figure 1, inconsistent with the pronounced "dip" of the CZ model. Using the LCSR approach it was also shown for many examples, see, e.g., refs. [17][18][19][20][21][22][23][24][25], that the CZ model leads to very large soft contributions to hard reactions, which contradict the data. Nevertheless, the paradigm "asymptotic-like LCDA versus CZ-like LCDA" continues to be the preferred language of many model studies.
A new wave of interest in the pion LCDA was inspired by the BaBar measurement [26] of the pion transition form factor γγ * → π that indicates very strong scaling violations up to the highest virtualities Q 2 ∼ 40 GeV 2 available. In order to explain this behavior, an unconventional "constant" shape of the pion LCDA was proposed [27,28], which triggered further discussion, see, e.g., ref. [22]. Although the similar Belle experiment [29] does not suggest strong scaling violations, the problem is far from being resolved and this measurement will be repeated by Belle II at the upgraded SuperKEKB accelerator at KEK [30] with a much improved projected precision. Motivated by these experimental needs and in the absence of a convincing first-principles calculation, the pion LCDA continues to attract a lot of attention. In the last few years several new calculations appeared, most notably using techniques based on Dyson-Schwinger equations (DSE) [31]. A short overview of several existing models and their distinctive features can be found in ref. [32]. For further models see, e.g., refs. [33,34].
Within the past 10-20 years lattice QCD has firmly established itself as the method of choice for nonperturbative calculations in QCD, as it has the potential to provide quantitative results with full control over all sources of uncertainty. The problem that we address here, however, is not simple. Lattice calculations of moments of the pion LCDA were proposed more than 30 years ago [35,36]. First pioneering studies were carried out within the quenched approximation in refs. [36][37][38][39] and with N f = 2 Wilson fermions in ref. [40]. The first modern calculations were performed more than a decade ago by the QCDSF/UKQCD collaboration using N f = 2 nonperturbatively improved Wilson fermions [41] and somewhat later by RBC/UKQCD [42] as part of their N f = 2 + 1 domain-wall fermion phenomenology program. More recently, the study of ref. [41] was extended in ref. [43] to a larger set of lattice ensembles with different volumes, lattice spacings, and pion masses down to m π = 150 MeV, also implementing several technical improvements. In this way the errors due to the chiral extrapolation could be brought under control but still no controlled continuum limit extrapolation could be carried out.
In this paper we close this last gap and present results of the first lattice calculation of the two lowest moments of the pion and kaon light-cone distribution amplitudes with full control of all systematic errors. This progress has become possible by the CLS (Coordinated Lattice Simulations) community effort [44] aiming at the production of very fine lattices using open boundary conditions in time and further algorithmic improvements to reduce the autocorrelations within the Monte-Carlo time-series. (Autocorrelations increase as the continuum limit is approached.) The calculation reported in this work has been carried out on 35 ensembles (see appendix A for details) using N f = 2 + 1 flavors of nonperturbatively improved Wilson (clover) fermions with pion masses down to the physical point, employing 5 different lattice spacings down to a = 0.039 fm. In addition, we use the momentum smearing technique [45], which enables us to reduce statistical fluctuations by improving the overlap of the meson interpolating field with the ground state. Employing this technique, first results for the second moment of the pion LCDA for a single lattice spacing were reported in ref. [46]. Since then we have enlarged the operator basis (cf. also ref. [47]) and added four lattice spacings as well as other quark mass combinations. The results are then obtained pursuing combined chiral and continuum limit extrapolations, utilizing data from three separate trajectories in the quark mass plane. As a by-product we also obtain the continuum limit quark mass dependence of the LCDA moments. A similar determination of the wave function normalization constants and the first LCDA moments of the lowest-lying baryon octet can be found in the companion article [48]. This article is organized as follows. In section 2 we first introduce LCDAs as well as the operators and correlators used in our analysis. Next, the renormalization of the lattice matrix elements is explained. This includes two steps: nonperturbative renormalization in the RI /SMOM (or RI /MOM) scheme and perturbative conversion from this scheme to the MS scheme. In section 3 we describe the set of gauge ensembles employed. Subsequently, we detail the analysis of the correlation functions (including the specific choice of operators and external momenta) and extract the relevant matrix elements from the lattice. We also provide the extrapolation formulae for the quark mass and lattice spacing dependence. Finally, in section 4, we present our results for the LCDA moments and assess the error budget, before we discuss our findings and confront these with values from the literature in section 5.

Continuum definitions
Each pseudoscalar meson has only one independent leading-twist LCDA, φ M , which can be defined via a meson-to-vacuum matrix element of a renormalized nonlocal quark-antiquark light-ray operator, where we consider the pion (M = π + ) withq =d and the kaon (M = K + ) withq =s. Here, z 1,2 are real numbers, n µ is an auxiliary light-like (n 2 = 0) vector, and |M (p) represents the ground state meson M with on-shell momentum p 2 = m 2 M . The light-like Wilson line connecting the quark fields, [z 2 n, z 1 n], is inserted to secure gauge invariance. The scale dependence of φ M is indicated by the argument µ 2 . We denote the quark masses as m q .
Neglecting the isospin breaking due to electromagnetic effects and nondegenerate light quark masses (by setting m u = m d ≡ m ), the LCDAs of all (charged and neutral) pions are trivially related such that it is sufficient to consider only one representative; the same holds for the kaons. The decay constant f M appearing in eq. (2.1) can be obtained as the matrix element of a local operator, and has the value f π ≈ 130 MeV [49] for the pion and f K ≈ 156 MeV [50] for the kaon. Within eq. (2.1) a fraction x of the longitudinal meson momentum is carried by the u quark, while theq antiquark carries the remaining fraction 1 − x. The difference of the momentum fractions is usually denoted as 3) The complete information on the LCDA can be encoded in a set of moments. One such set is defined by Another possible set of moments is a M n (µ 2 ) = 2(2n + 3) 3(n + 1)(n + 2) where C 3/2 n (ξ) are Gegenbauer polynomials, which correspond to irreducible representations of the collinear conformal group SL(2, R). Both sets, the ξ-moments ξ n and the Gegenbauer moments a M n , are related by a simple linear transformation, cf. eqs. (2.15b) and (2.16b) below. 1 Since the Gegenbauer polynomials form a complete set of functions, the LCDAs can be expanded as where the coefficients a M n are renormalized multiplicatively in leading logarithmic order as a consequence of conformal symmetry [51]. Due to C-parity, all odd moments of the pion, i.e., ξ n π and a π n for n = 1, 3, . . . , vanish in the limit of exact isospin symmetry. Higherorder contributions in the Gegenbauer expansion are suppressed at large scales, since the anomalous dimensions of a M n increase with n [5]. Hence, in the asymptotic limit µ → ∞ only the leading term survives, which is usually referred to as the asymptotic LCDA. From here on we will suppress the explicit scale dependence of the DAs and their moments in the notation. Our lattice results will be given at the fixed scale µ = 2 GeV in the MS scheme with three active flavors.

Lattice definitions
From now on we will work in Euclidean spacetime and follow the conventions of ref. [43]. The renormalized light-ray operator on the left-hand side of eq. (2.1) generates renormalized local operators. This means that the moments (2.4) of the LCDAs can be expressed in terms of matrix elements of local operators that can be evaluated using lattice QCD. In order to calculate the first and second moments of the pseudoscalar LCDAs we define the bare lattice operators where the covariant derivative D µ is discretized symmetrically. To obtain a leading-twist projection we symmetrize over all Lorentz indices and subtract all traces. This procedure is indicated by parentheses, e.g., In principle, one could also consider an operator O + ρµ , replacing the minus sign in eq. (2.8c) by a plus sign. However, as O + ρµ differs in C-parity from O − ρµ , these two operators cannot mix with each other so that O + ρµ is irrelevant for our calculation. In contrast, the operator O + ρµν has the same C-parity as O − ρµν and must be taken into account. Introducing the shorthand notation ⃡ On a hypercubic lattice, the continuous O(4) symmetry is reduced to its discrete H(4) subgroup. This symmetry breaking can in principle induce mixing of the operators of interest with lower-dimensional operators accompanied by coefficient functions that diverge with a power of 1/a. For the first two ξ-moments this mixing can be avoided by selecting lattice operators that belong to a suitable irreducible representation of the hypercubic group H(4) [41,42]. For the calculation of the first moment we use the operators O − 4µ , while for the second moments we choose O ± ρµν with all three indices different, see also section 2.3.
In order to extract the desired moments we use two-point correlation functions of the operators with an interpolating current, where the index p indicates that the quarks appearing within the interpolator (2.8a) have been momentum smeared [45,46] (employing APE smeared [52] spatial gauge transporters) to optimize the overlap with the ground state. The smearing parameters are not only adjusted according to the momentum but also optimized with respect to lattice spacing and quark mass. The ground state will dominate for sufficiently large values of the sourcesink separation t. In this limit, neglecting effects from the temporal boundaries, one obtains For the extraction of the first moment we consider the ratios Similarly, for the required matrix elements for the second moment we consider In contrast to the ratios (2.11), the two ratios defined in eqs. (2.12) transform according to the same irreducible representation of H(4) and will give the same result R ± 2 = R ± 2,a 1 = R ± 2,a 2 (in the limit t → ∞, p j a −1 ). However, R ± 2,a 1 and R ± 2,a 2 are affected differently by excited states, cf. section 3.2.

Renormalization procedure
The lattice operators have to be renormalized to obtain matrix elements in the MS scheme. As mentioned above, the continuous Euclidean O(4) symmetry is reduced to that of its finite hypercubic subgroup H(4) on the lattice. Therefore, symmetry imposes much weaker constraints on the mixing of operators under renormalization. In order to avoid mixing as far as possible, in particular mixing with lower-dimensional operators, we use operators from suitably chosen multiplets that transform according to irreducible representations of H(4) and possess a definite C-parity. In the case of the operators (2.8c) with one derivative we consider two multiplets transforming according to nonequivalent representations: one, labeled a, consisting of the six operators O − ρµ with 1 ≤ µ < ρ ≤ 4 and another one, The operators (2.8d) and (2.8e) with two derivatives have equal C-parity and behave identically under both continuum and lattice spacetime transformations. Hence, they will necessarily mix with each other. We utilize the multiplets which transform under H(4) according to one and the same four-dimensional irreducible representation. Their symmetry properties guarantee that they do not mix with any other operators. We determine the renormalization and mixing coefficients nonperturbatively on the lattice using the same RI /SMOM scheme [53] as was used in ref. [43]. For the coarser lattice spacings (β = 3.4, 3.46, 3.55) we have ensembles with different quark mass values m = m s and (anti-)periodic boundary conditions in time at our disposal so that we can proceed in exactly the same way as in ref. [43], starting from Landau-gauge-fixed threepoint functions where O represents the operators from eqs. (2.8b)-(2.8e) with an antiquark flavorq =d that is mass-degenerate with the u quark. However, a problem arises on the finer lattices. For β = 3.7 and 3.85 we are forced to work with open boundary conditions in time to reduce autocorrelations in the Monte-Carlo time-series [54,55]. In this case we modify the computation of the required three-point functions in two respects: we place the momentum sources within a subvolume, keeping a sufficiently large distance from the boundaries in the time direction, and we restrict the (final) sum over z to an even smaller volume inside this subvolume. The further analysis can then be performed as in the periodic case. A detailed discussion, including a justification of this method and a comparison with the results from periodic boundary conditions, will be the topic of a dedicated, forthcoming publication. The ensembles with symmetric quark masses (m = m s ) used for the calculation of the renormalization factors are detailed in table 8. Unfortunately, we could only afford to generate ensembles for two distinct values of m = m s at β = 3.7 and 3.85. In the other cases the mass dependence of the amputated three-point functions is rather mild, so that we are confident that this restriction does not significantly affect the reliability of the required chiral extrapolations.
In the case of the first LCDA moment of the kaon it is also possible to carry out the renormalization via the RI /MOM scheme [56,57], where even the three-loop matching to the MS scheme is available [58][59][60]. Therefore, we choose to present four distinct results: with one-and two-loop matching [61,62] via the RI /SMOM scheme as well as with twoand three-loop matching using the RI /MOM scheme.
The tiny statistical errors of the results are negligible in comparison to the systematic uncertainties. In order to estimate the latter we proceed similarly to ref. [43] and perform a number of fits, varying one element of the analysis at a time. We carry out two independent determinations of the renormalization and mixing coefficients, namely with one-loop and two-loop truncations of the perturbative expansion of the conversion factors from the RI /SMOM scheme to the MS scheme for use in NLO and NNLO calculations in perturbative QCD, respectively. In both cases we vary the initial scale µ 1 of the fit range and the number n disc of terms in the parametrization A 1 a 2 µ 2 + · · · + A n disc (a 2 µ 2 ) n disc of the lattice artifacts. In order to take into account the uncertainties in the determination of the lattice spacing, the central values of 1/a 2 shown in table 1 are multiplied by a factor λ 2 scale = 1.03. This value contains the scale uncertainty of 8t * 0 = µ * −2 ref given in ref. [63] and the largest error of our determination of t * 0 /a 2 , added in quadrature. Finally, also Λ We determine the LCDA moments separately for each of the resulting renormalization and mixing coefficients, thus generating a set of five values per renormalization scheme at a given loop order. In this way we obtain two sets of results for the second LCDA moments, one using the two-loop SMOM conversion factors and another one employing the one-loop SMOM conversion factors. As explained above, for the first moment of the kaon LCDA we even have four such sets of results, as we can also nonperturbatively convert the bare lattice results to the RI /MOM scheme instead and then utilize the two-loop or three-loop conversion factors between the RI /MOM and the MS schemes.
In each set we take the results of fit 1 as our central values. Defining δ i , i = 2, 3, 4, 5, as the difference between the number based on fit i and the result based on fit 1, we estimate the systematic uncertainties due to the renormalization factors as The dominant uncertainties are related to the low-momentum cut-off of our fit range (δ 2 ), i.e., the scale dependence, and the parametrization of lattice artifacts (δ 3 ). The former becomes smaller when going from one-loop to two-loop perturbative accuracy, while the latter uncertainty shrinks as the lattice spacing is reduced. The uncertainty induced by the scale setting (δ 4 ) and the error of the strong coupling parameter (δ 5 ) are negligible. Note that all figures in this article showing renormalized data are generated using the RI /SMOM intermediate scheme with two-loop matching to the MS scheme.
Finally, the renormalized first moments are related to the ratios defined in eqs. (2.11) by while the second moments are related to the ratios (2.12) via In the continuum 1 2 MS = 1, while it can differ from unity on the lattice, see section 4.1.
The ζs denote ratios of the renormalization constants of the operators (2.8c)-(2.8e) over the renormalization constant of the axialvector current (2.8b), cf. ref. [43]. Henceforth, ξ n , 1 n , and a n are always implied to be renormalized in the MS scheme and we omit the superscript MS.  Figure 2. Schematic illustration of the mass trajectories of the lattice ensembles used in this study. Along the flavor symmetric line (blue) all pseudoscalar mesons have equal mass (m 2 K = m 2 π ), which is equivalent to equal quark masses (m = m s ). The (green) line of the physical average quadratic meson mass (2m 2 K + m 2 π = phys.) corresponds to an approximately physical mean quark mass (2m + m s ≈ phys.). The red line is defined by 2m 2 K − m 2 π = phys. and indicates an approximately physical strange quark mass (m s ≈ phys.). The gray dot marks the physical point.

Lattice ensembles
We use lattice ensembles generated within the CLS effort [44] employing N f = 2 + 1 flavors of nonperturbatively O(a) improved Wilson fermions [64,65] combined with the tree-level Symanzik improved gauge action [66]. For details on the action and the simulation see ref. [44]. 2 Since that publication more CLS simulation points have been added, see, e.g., ref. [68]. An overview of the ensembles analyzed here is given in appendix A. Most CLS ensembles use open boundary conditions in the time direction, which allows us to carry out simulations at very fine lattice spacings without facing the problem of topological charge freezing [54,55].
Five values of the inverse coupling constant β = 6/g 2 are realized, corresponding to lattice spacings ranging from a = 0.086 fm down to a = 0.039 fm, see table 1. Here we set the scale using √ 8t * 0 = 0.413(6) fm [63], where t * 0 is defined in ref. [69] as the Wilson flow scale t 0 [70], computed at a particular reference point in the quark mass plane. The numerical value was obtained by matching the average continuum limit pion and kaon decay constant f πK = (2f K + f π )/3 to experiment [69].
At each lattice spacing we have several points in the quark mass plane, along three trajectories: (a) along a nearly-physical fixed value of the trace of the mass matrix Tr M ≡ m u + m d + m s = 2m + m s = phys., (b) varying the light quark mass while trying to keep the renormalized strange quark mass m s constant at its physical value, and (c) along the "symmetric" line m = m s , where light and strange quark masses are equal. The first two trajectories intersect close to the physical quark mass point. The locations of these three lines are shown in figure 2. We determine the LCDA moments on various ensembles along these trajectories; our largest pion mass is about 420 MeV and the smallest one is 130 MeV. Table 6 of appendix A contains all lattices lying on line (a) (Tr M = constant). This line starts with a lattice at the flavor symmetric point and approaches the physical point, decreasing the light quark mass while simultaneously increasing the strange quark mass.  The spatial extents of the lattices used to determine the LCDA moments are always larger than 2.4 fm and, with very few exceptions, larger than four times the inverse mass of the lightest pseudoscalar meson, see also tables 6-8. For the pseudoscalar meson masses the expected corrections due to finite volume effects calculated at next-to-leading order in chiral perturbation theory (ChPT) [71,72] are smaller than half of their statistical errors. To this order the LCDAs are not affected by finite volume corrections at all since they are normalized with respect to the decay constant, see eq. (2.1). Therefore, it is well justified to neglect volume effects in our analysis.

Analysis of correlation functions
Below we specify our choice of correlators and momentum directions. For the first moment we have operators from two different H(4) multiplets at our disposal (cf. eqs. (2.11)). For the ratio in eq. (2.11a) we select the momenta p = (±1, 0, 0)℘, p = (0, ±1, 0)℘, and p = (0, 0, ±1)℘, where ℘ = 2π L . We then extract R − 1,a as a function of t according to where the forward/backward momentum averaging is performed by the operatorp ± : For the ratio in eq. (2.11b) we may simply set p = 0 to obtain We then renormalize the above ratios, multiplying by ζ a and ζ b according to eq. (2.15a). Finally, ξ 1 K is obtained by carrying out a simultaneous fit to the plateau that is reached at large t-values as depicted in figure 3.
For the extraction of the second moments one needs at least two nonvanishing momentum components, cf. eqs. (2.12). We have already addressed the problem of the deterioration of the signal-to-noise ratio for increasing momenta |p| in our previous work [46], where we proposed to employ the momentum smearing technique (introduced in ref. [45]) for all quark sources in order to improve the statistical error and to reduce contributions from excited states. The momentum smearing technique requires two inversions per momentum direction and in order to evaluate the full sum in eq. (2.12a) we performed six inversions to realize the momenta p = (1, 1, 0)℘, p = (1, 0, 1)℘, and p = (0, 1, 1)℘ in ref. [46]. In the present work we select the slightly higher momentum p = (1, 1, 1)℘, which allows us to evaluate both eqs. (2.12a) and (2.12b). This requires only two inversions in total. We compare the two ratios R ± 2,a 1 and R ± 2,a 2 for this momentum in figure 4. We see that R + 2,a 2 is by far superior for the extraction of R + 2 , while R − 2,a 1 is preferable for the determination of R − 2 . Since the operators O 4ij and O 123 belong to the same H(4) multiplet, combining the results for R + 2,a 2 and R − 2,a 1 in order to obtain ξ 2 via eq. (2.16a) is allowed and does not require any additional considerations regarding the renormalization.
As shown in [46], larger momenta can even improve the signal-to-noise ratio in certain situations. This is not the case here: the correlation functions with p = (1, 1, 1)℘ have a slightly inferior signal-to-noise ratio compared to those using p = (1, 1, 0)℘, cf. eq. (27) of ref. [46]. However, this choice enables us to obtain results for the whole operator multiplets in eqs. (2.13) from a single momentum, which makes the calculation more efficient (roughly by a factor of four). That the additional ratio R + 2,a 2 yields a much better ground state plateau (see the left panel of figure 4) is an extra benefit.   = (1, 1, 1)℘ for the ensemble N203 (a = 0.064 fm). Clearly, for the extraction of R + 2 (left), the ratio R + 2,a2 is to be preferred to R + 2,a1 , which suffers considerably from excited state effects and carries larger statistical errors. For the case of R − 2 (right) neither data set seems to indicate any significant excited state contribution, but the statistical errors of R − 2,a1 are much smaller. The bands indicate the fit ranges and results.

Chiral extrapolation
The CLS ensembles described in section 3.1 (for more detail see appendix A) enable us to perform a joint chiral and continuum limit extrapolation. As will be explained in section 4, both limits are well controlled, the latter due to the extended set of different lattice spacings at our disposal and the former due to the approach of the physical point along two distinct quark mass trajectories, with further constraints from the points along the symmetric line. The formulae for the chiral extrapolation of the first two LCDA ξ-moments of the lowestlying pseudoscalar meson octet, i.e., the π, the K, and the η 8 mesons, 3 have been worked out in ref. [73]. For the even moments ξ 2n M one obtains where α (2n) and β (2n) are low energy constants (LECs). It is convenient to introduce the variables such thatm 2 is approximately constant along the Tr M = constant trajectory, while δm 2 vanishes for degenerate quark masses m = m s . Here B 0 = | ūu |/F 2 0 ≈ 2| ūu |/f 2 π is the quark condensate parameter. Along the symmetric line the mesons have to form an exact SU(3) flavor octet with one and the same leading-twist LCDA for the π, the K and the η 8 . This becomes evident when rewriting eqs. (3.4) in terms of the new variables: Here,Ā (2n) = 2α (2n) + 3β (2n) /(2B 0 ) and δA (2n) = α (2n) /(3B 0 ) are linear combinations of the LECs of eqs. (3.4). Note that the breaking of SU(3) flavor symmetry is highly constrained as, to one-loop order in ChPT, we have only one independent symmetry breaking parameter δA (2n) per LCDA moment. This will allow us to infer the shape of the η 8 LCDA from the pion and kaon data.
In the limit of exact isospin symmetry, C-parity implies that the LCDAs of the pion and η 8 are even functions of ξ. Therefore, the odd moments vanish. This also applies to the LCDA of the kaon in the limit of exact flavor symmetry δm 2 = 0. Therefore, re-expressing the corresponding formulae of ref. [73] in terms of the variablesm and δm gives for the odd moments

Discretization effects
For both LCDA moments we expect the leading-order discretization effects to be linear in a, as the corresponding operators, O − ρµ and O − ρµν , have not been O(a) improved. 4 We make the ansatz where the chiral extrapolation formulae of section 3.3 are combined with a linear parametrization of discretization effects, including mass-dependent terms. The SU(3) flavor constraints will be violated by O(a) terms since our fermion formulation explicitly breaks chiral symmetry. Therefore, δc (2) π and δc (2) K are independent parameters. Within this ansatz we require a total of four parameters to describe the lattice spacing and quark mass dependence of ξ 1 K , while seven parameters are needed for our joint extrapolation of ξ 2 π and ξ 2 K that also yields ξ 2 η 8 . We will see that all lattice data are well described by the above ansätze. Nevertheless, we will vary the parametrization to explore the systematics associated with the choice of this particular functional dependence.
In the continuum, the remaining operator O + ρµν can be written as the second derivative of the axialvector current, . This is not the case on the lattice and the renormalization factors of O + and A differ. However, in the continuum limit the renormalized lattice ratio should approach unity, such that the continuum relation a 2 = 7 12 5 ξ 2 − 1 is recovered from eq. (2.16b). We employ a nonperturbatively O(a) improved fermion action and tree-level O(a) improved derivatives in our operators. Assuming small order a discretization effects in 1 2 M , 0,2 a 2 +ē (2) 2m 2 a 2 + δe (2) M,2 δm 2 a 2 (3.10) should provide a sensible parametrization of the data. In the next section we will discuss and check this ansatz.

Extrapolation strategy and error budget
In the following we present our results for the first and the second ξ-moments and Gegenbauer moments of the leading-twist pseudoscalar meson distribution amplitudes. In addition to the results for the pion and the kaon, which are extracted directly from the lattice data, we infer the second moment of the η 8 meson using eq. (3.6c) from the SU ( Figure 5. The quantity 1 2 M as a function of the squared lattice spacing a 2 , plotted at the physical mass point. The solid lines represent the result of a global fit using eq. (4.1). The points shown have been obtained by translating all data points to the physical masses along the fitted function and then averaging measurements from the same lattice spacing. The dashed curves correspond to the alternative fit carried out to investigate linear terms as described in the main text.
breaking constraints obtained from ChPT in ref. [73]. Previous lattice determinations of the Gegenbauer moments [36-43, 46, 74-78] lacked ensembles with lattice spacings smaller than 0.06 fm and so far no controlled continuum limit extrapolation has been carried out. This is particularly problematic for the second moment a M 2 , which mixes with 1 2 M under renormalization, see eqs. (2.16).

A game of ones
The continuum limit 1 2 M a→0 − −− → 1 is known. Using this value as a constraint and fitting our data, we find that the a dependence is mostly quadratic and the possible linear contribution is small. This is consistent with expectations based on tree-level lattice perturbation theory, where linear terms vanish exactly.
One can play another game, pretend that the continuum value of 1 2 M is not known, and try to determine it from the data. The quadratic fit ansatz 0,2 a 2 +ē using I M as a free parameter, gives a continuum limit value close to one with only 0.5% deviation, see the solid line in figure 5. This agreement is nontrivial (unrenormalized lattice values in the considered region of lattice spacings lie in the range 0.59-0.68, see, e.g., the left panel of figure 4) and can be viewed as confirmation of our calculation of the corresponding renormalization constant. However, without the constraint at a = 0, the smallness of linear contributions in comparison to the quadratic a dependence cannot be inferred from the data: an alternative fit including the additional linear terms e (2) 0 a,ē (2)m2 a, and δe (2) M δm 2 a (dashed curve in figure 5) leads to a continuum value that is about 2.5% above unity. The difference can be viewed as a systematic uncertainty of the continuum extrapolation (labeled a in the following), yielding the "lattice values" I π = 0.9947 +2 −2 (80) r (301) a and I K = 0.9941 +1 −2 (80) r (300) a , where statistical errors are given by the sub-/superscript pair and the uncertainty due to the renormalization (r) is determined as described in section 2.3. To avoid misunderstanding: the values of I M (and the fits shown in figure 5) are not used in the determination of the moments of meson LCDAs, to be discussed in the following sections. Their determination merely serves as a sanity check to strengthen the confidence in our renormalization procedure.  Figure 6. Dependence of the moments ξ 2 M on the squared pion mass, plotted in the continuum limit. The points shown have been obtained by translating all data along the fitted function (keeping the masses fixed). The plots for the individual lattice spacings can be found in figure 12 in appendix A. The solid lines and shaded statistical error bands represent our main result. The dashed curves correspond to the mean value of an alternative fit (including a term of higher order in the masses) used to estimate the parametrization dependence as described in the main text.
In comparison to our previous work, see figure 3 of ref. [43], we achieve a much higher statistical precision for 1 2 M , such that the statistical error now contributes by far the smallest uncertainty. This improvement in statistics is mostly due to employing the operator O + 123 in the new method (2.12b), compared to the old method involving the operators O + 4ij , see also the left panel of figure 4. Furthermore, it turns out that also the systematic uncertainties due to renormalization (0.8%) and due to discretization effects (3%) are quite small.

Extrapolation of the second LCDA moments
For the extrapolation of the second moments ξ 2 π and ξ 2 K we use eq. (3.8b). We then insert the fitted LECs ξ 2 0 ,Ā (2) , and δA (2) into eq. (3.6c) in order to obtain a prediction for ξ 2 η 8 in the continuum. The combined extrapolation is shown in figure 6 as a function of the pion mass and in figure 7 as a function of the lattice spacing. Figure 6 shows that the breaking of SU(3) flavor symmetry among these observables is rather small. Actually, within our errors, we find no differences between ξ 2 π , ξ 2 K , and ξ 2 η 8 . To estimate the systematic uncertainty due to the choice of the parametrization of the mass dependence we perform an alternative fit by including the additional termĀ  figure 6 and we take the difference with respect to the mean value of our main fit as the corresponding error. It can be seen that the second moments of the pseudoscalar LCDAs depend only mildly on the quark masses. In contrast, the discretization effects are quite sizable and amount to a correction of roughly 10% from our largest lattice spacing of a = 0.086 fm to the continuum, as shown in figure 7. To estimate the systematics of the a dependence we again perform an alternative fit, this time adding the term c    Figure 7. Dependence of the moments ξ 2 M on the lattice spacing a, plotted at physical quark masses. The points shown have been obtained by translating all data along the fitted function (keeping the lattice spacing fixed) and then averaging measurements with the same a. The plots for the individual trajectories can be found in figure 12 in appendix A. The solid lines and shaded statistical error bands represent our main result. The dashed curves correspond to the mean value of an alternative fit (including a term proportional to a 2 ) used to estimate the parametrization dependence as described in the main text.
We have checked that other methods to estimate this systematic error lead to compatible results, e.g., omitting the data from the coarsest lattice spacing. Another possibility is to consider continuum extrapolations for two lattice observables that have the same continuum limit. To this end, we compare the second Gegenbauer moment a M 2 defined in terms of ξ 2 M via the continuum theory relation which is natural at a finite lattice spacing. As argued in section 4.1, the difference between these two quantities should be mainly due to O(a 2 ) effects. A comparison is shown in figure 8 for the pion (left) and kaon (right). In both cases we perform a linear extrapolation in the lattice spacing. The difference in the continuum compares reasonably well to the estimates for O(a 2 ) effects obtained from the procedure explained above.

Extrapolation of the first LCDA moment
A combined continuum and chiral extrapolation of ξ 1 K is performed using eq. (3.8a), which automatically enforces the constraint that all odd moments have to vanish in the limit of exact flavor symmetry (which is also true for the lattice data). We therefore only have data points for lattices with nondegenerate quark masses, see figure 9 (left). The mass dependence in the continuum limit is determined by the single parameter δA (1) = 0.141(23) GeV −2 , see eq. (3.8a). Notably, we find only a very mild dependence of the first moment on the lattice spacing that is consistent with a flat behavior within errors, see figure 9 (right). The parametrization dependence is investigated, as above, by performing two alternative fits, each including a single additional parameter. These fits are indicated  Figure 9. Left: The same as figure 6, but for the moment ξ 1 K ; the two relevant trajectories have been condensed into one plot. Right: The same as figure 7, but for the moment ξ 1 K . The plots for the individual lattice spacings and trajectories can be found in figure 13 in appendix A. by dashed lines in the corresponding plots; one includes the term δĀ (1)m2 δm 2 for the mass dependence, 6 the other one includes the term c (1) 0,2 a 2 for the lattice spacing dependence.

Summary of the results
The mass dependence of the first two Gegenbauer moments a M 1 = 5 3 ξ 1 M and a M 2 = 7 12 5 ξ 2 M − 1 in the continuum limit is summarized in figure 10. Our final results for 6 One could use a term ∝ δm 4 instead (this adds a single parameter in the case of the odd moments), which leads to a very similar estimate for the uncertainty. Usingm 4 is however not allowed since the whole fit function must be proportional to δm 2 due to symmetry.
the moments ξ 1 M and ξ 2 M as well as the corresponding Gegenbauer moments (in the continuum limit at the MS scale µ = 2 GeV) are collected in table 3. It can be seen as a success of our strategy, i.e., generating ensembles on different quark mass trajectories while simultaneously reaching fine lattice spacings, that all the systematic uncertainties can be controlled and are of a similar or smaller size than the statistical accuracy. In analogy to the prevalent procedure used in determinations of parton distribution functions from experimental data, we quote separate results for the NLO (one-loop) and the NNLO (two-loop) analysis. Even though the results obtained using the SMOM scheme with NLO and NNLO matching almost agree within the given renormalization error, the central values still deviate considerably from each other so that a three-loop matching formula between the RI /SMOM and MS schemes would be welcome. As to be expected, the systematic uncertainty due to renormalization decreases for increasing loop order. We quote our SMOM NNLO values as the final results in the abstract.

Discussion
In table 4 we compare our result for the second moment of the pion LCDA to values from the literature. Our number is compatible with the previous result [43] obtained several years ago with N f = 2 clover fermions. 7 The quality of the present data is much higher, enabling a controlled continuum extrapolation with quantifiable errors. Our result for a π 2 is smaller by a factor of four in comparison to the original CZ calculation [8,10] evolved to the 2 GeV scale, but the difference to more recent QCD sum rule calculations is much smaller and in particular the sum rules involving nonlocal vacuum condensates [14,15,25,83] yield an estimate that is consistent with our results within the quoted error bar. The entries in table 4 marked "LCSR" are obtained from experimental data in the factorization framework using LCSR-corrected coefficient functions to take into account the contributions of "soft" regions. It is interesting that new data from the BaBar [26] and Belle collaborations [29] generally support small values of the second moment, compatible with our result. Methods based on Dyson-Schwinger equations (DSE) [31] suggest somewhat larger values.
A similar comparison for the first two moments of the kaon is presented in table 5. Our result for the first moment is consistent with earlier lattice calculations as well as with results from QCD sum rules and is somewhat smaller compared to the DSE calculation in ref. [94]. Regarding the second moment of the kaon LCDA, our number is lower than "old" lattice estimates [41,42,77] but agrees remarkably well with the DSE prediction [94] based on the so-called DCSB-improved version of the truncation.
As far as future calculations of the second moment of the pion and kaon LCDAs are concerned, the accuracy can be improved by increasing the statistics in particular for the ensembles at small lattice spacings and quark masses but also by adding additional simulation points. Also a three-loop calculation of the perturbative matching to the RI /SMOM scheme is required to improve the overall accuracy.
Regarding phenomenological applications, the first inverse moment which is equal to the sum of all Gegenbauer coefficients, is of particular importance since this quantity enters at leading order in factorization theorems (see, e.g., ref. [10]). Unfor- 7 The result of ref. [43] does not correspond to the continuum limit but to an average of data within a window of lattice spacings a ≈ 0.06-0.08 fm. Moreover, in this reference the values of a2 and ξ 2 are related via eq. (4.3), where 1 2 = 1 for a > 0. Directly comparing our results to those of ref. [43] at a finite lattice spacing may be misleading as in that simulation a different number of sea quarks and a different gluonic action were used. Table 4. The second moment of the pion LCDA at the MS scale µ = 2 GeV. The CZ model fixes a π 2 = 2/3 at the low scale µ 500 MeV; for a discussion of the extrapolation to higher scales see ref. [79]. The abbreviations stand for: LQCD: lattice calculation; N f = 2(+1): calculation using N f = 2(+1) sea quarks; SW: nonperturbatively O(a) improved Sheikholeslami-Wohlert (i.e., Wilson-clover) fermion action; DWF: domain-wall fermions; QCDSR: QCD sum rules; NLC: nonlocal condensates; LCSR: light-cone sum rules; R: renormalon model for twist-4 corrections; DSE: Dyson-Schwinger equations with rainbow-ladder truncated (RL) or DCSB-improved (DB) kernels. The LCSR analysis is based on the experimental data from the CLEO [80], BaBar [26], and Belle [29] collaborations. Among previous lattice studies only in ref. [41] an attempt of a continuum limit extrapolation was made. The result of ref. [39]  tunately, there is no known way to evaluate it directly on the lattice. As an illustration, we compare two phenomenologically acceptable models of the pion LCDA. The first model is the expansion in Gegenbauer polynomials truncated after n = 2 and the second model is based on a simple power-law parametrization: φ as (x) φ as (x) Figure 11. The truncated Gegenbauer expansion (5.2a) and the power-law parametrization (5.2b) at µ = 2 GeV obtained using our results for a M 1 and a M 2 . The left panel shows the resulting DAs for the pion, which are symmetric under x ↔ 1 − x, while the results for the kaon DA on the right panel are slightly skewed towards the strange quark due to flavor symmetry breaking. In all cases the deviation from the asymptotic shape is significant.
Both formulae have two parameters, where for the pion, of course, a π 1 = 0 and α − π = α + π . We fix them such that our calculated values for a M 1 and a M 2 (in the SMOM scheme at two-loop order, cf. table 3) are exactly reproduced also in the second model. Hence, both models have by construction the same value for the first two Gegenbauer coefficients, but differ in higher-order coefficients.
The results are shown in figure 11. Both models are somewhat "flatter" in comparison to the asymptotic LCDA shown by the gray curve, and in general do not seem to differ very much. The model dependence of the first inverse moment is, however, sizable. We obtain for the pion where the errors have been obtained by adding the individual errors of table 3 in quadrature.
Both numbers are phenomenologically viable, in particular the second one is very close to (1 − x) −1 π = 3.6 (at µ = 2 GeV) from the model of ref. [86], which provides a good description of the Belle data [29] for the πγγ * form factor.
The QCD description of form factors based on our models (I) and (II) will differ by as much as 10%. The necessity to go beyond the second Gegenbauer moment is thus obvious. A brute-force extension of the present approach to operators with a larger number of derivatives does not seem to be viable even if the problem of the mixing with lowerdimensional operators is solved. Consider in particular the fourth moment, ξ 4 M = 3/35 + (8/35)a M 2 + (8/77)a M 4 , for which we obtain in the two models One sees that even if both ξ 2 π and ξ 4 π were measured with 1% precision on the lattice (which is already optimistic given our statistical error of ∼ 2.5% on ξ 2 π ), the value of a π 4 cannot be extracted reliably as it is overshadowed by the uncertainty in a π 2 . Therefore, alternative methods should also be investigated.
In the past few years exploratory studies appeared aiming at the extraction of the pion LCDA from lattice calculations of suitable Euclidean correlation functions in position space [95][96][97][98], see also refs. [99][100][101]. After taking the continuum and other appropriate limits, these can be expressed in terms of LCDAs in the framework of QCD factorization within the continuum theory, in analogy to the extraction of parton distributions from fits to experimentally measured structure functions. In other words, the role of lattice QCD is in this case to provide a complementary set of observables from which the LCDAs can be extracted. In particular, in ref. [98] it has been demonstrated that using the approach of ref. [95], the contributions of different Gegenbauer moments can be separated, at least in principle, by considering the correlation functions at large "Ioffe times". These new techniques generally require hadron sources with very large momentum combined with good statistical accuracy and very fine lattices to control the corresponding discretization errors. Whether these position space methods or the moment method employed here will be more useful to constrain higher moments of LCDAs is at present unclear.
(grant nos. GA67-12, GA69-20, GA71-26), the PLGRID consortium for a computer time allocation on the Prometheus machine hosted by Cyfronet Krakow (grants hadronspectrum, nspt, pionda), PRACE (Partnership for Advanced Computing in Europe, http://www. prace-ri.eu) for awarding us access to the Marconi-KNL machine hosted by CINECA at Bologna, Italy, and the Leibniz Supercomputer Centre (LRZ, https://www.lrz.de) in Garching for access to the coolMUC3 cluster. Additional computations have been carried out on the Regensburg QPACE 2 computer [103] and the QPACE 3 machine of SFB/TRR-55. Some of the m = m s gauge ensembles used were generated by members of the Mainz group on the Wilson and Clover HPC Clusters of IKP Mainz. We acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time for GCS large-scale projects on the GCS share of the supercomputers SuperMUC at LRZ and JUQUEEN at JSC, where many of the ensembles used here were generated. GCS is the alliance of the three national supercomputing centers HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich) and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF).

A Lattice ensembles and supplementary figures
Below we list the properties of the analyzed lattice ensembles for the three quark mass trajectories: Tr M = phys. in table 6, m s = phys. in table 7, and m = m s in table 8. The latter also contains the ensembles that have been used solely for the determination of renormalization factors. We also show the results of the global fit for the second moments ξ 2 M in figure 12 and for the first moments ξ 1 M in figure 13. These are exactly the same fits that have been used to produce the more concise figures 6, 7, 9, and 10. In contrast to the figures of the main text we resolve the dependence on all relevant variables simultaneously, i.e., we display the full mass dependence along the three individual trajectories for each of the five lattice spacings as well as in the continuum limit. Table 6. List of the ensembles on the Tr M = phys. trajectory. The inverse gauge coupling β determines the lattice spacing (cf. table 1), while the spatial and temporal extents fix the lattice geometry N 3 s × N t . Boundary conditions in time direction are either periodic (p) or open (o). The light and strange hopping parameters, κ and κ s , determine the corresponding quark masses; the resulting approximate meson masses m π and m K are given in units of MeV, followed by the spatial lattice size in pion mass units. Finally, we give the number of gauge configurations used to measure the second moments.

Ens.
β N s N t bc κ κ s m π m K m π L conf.     Figure 13. The pion mass dependence of the moments ξ 1 M , defined in eq. (2.15a), plotted for all lattice spacings as well as in the continuum limit (where, for illustrative purposes, all points have been translated along the fitted function). The columns correspond to the three quark mass trajectories, cf. figure 2. The dotted gray lines mark the physical meson masses. Due to symmetry, this moment vanishes exactly for the π and η 8 mesons as well as on the symmetric line.