The pion and kaon $\langle x^3 \rangle$ from lattice QCD and PDF reconstruction from Mellin moments

We present a calculation of the pion and kaon Mellin moment $\langle x^3 \rangle$ extracted directly in lattice QCD using a three-derivative local operator. We use one ensemble of gauge configurations with two degenerate light, a strange and a charm quark ($N_f=2+1+1$) of maximally twisted mass fermions with clover improvement. The ensemble reproduces a pion mass $\sim260$ MeV, and a kaon mass $\sim530$ MeV. Excited-states contamination is evaluated using four values of the source-sink time separation within the range of $1.12-1.67$ fm. We use an operator that is free of mixing, and apply a multiplicative renormalization function calculated non-perturbatively. Our results are converted to the $\overline{\rm MS}$ scheme and evolved at a scale of 2 GeV, using three-loop expressions in perturbation theory. The final values are $\langle x^3 \rangle_\pi^{u^+}=0.024(18)_{\rm stat}(2)_{\rm syst}$, $\langle x^3 \rangle_K^{u^+}=0.035(6)_{\rm stat}(3)_{\rm syst}$, and $\langle x^3 \rangle_K^{s^+}=0.075(5)_{\rm stat}(1)_{\rm syst}$, where the systematic error is the uncertainty due to excited state contamination. We combine $\langle x^3 \rangle$ with the two lower moments to obtain the ratios $\langle x^3 \rangle/\langle x \rangle$ and $\langle x^3 \rangle/\langle x^2 \rangle$, as well as ratios between the pion and kaon moments. In addition, we reconstruct the $x$-dependence of the pion and kaon PDFs via 2- and 3-parameter fits to our results. We find that the reconstruction is feasible and that our lattice data favor a large $x$-dependence that falls as $(1-x)^2$ for both the pion and kaon PDFs. We integrate the reconstructed PDFs to extract the higher moments with $4\leq n\leq 6$. Finally, we compare the pion and kaon PDFs, as well as the ratios of their moments, to address the effect of SU(3) flavor symmetry breaking.


I. INTRODUCTION
The pion, kaon, and eta mesons comprise the octet of Nambu-Goldstone bosons, which are unique among hadrons because their masses vanish in the chiral limit. The valence quark structure of these mesons is given by a combination of a quark and an anti-quark with flavors u, d, and s, and if the masses of these quarks are equal then these Nambu-Goldstone bosons are mass degenerate (up to electroweak effects). Since the mass of the strange quark is significantly higher than that of the light u and d quarks -2 m s /(m u + m d ) = 27.46 ± 0.15 ± 0.41 [1] -comparison between pion and kaon observables provides a unique window into the interplay between strong interaction forces described by quantum chromodynamics (QCD) and quark mass effects [2].
For the pion and kaon, signficant SU(3) flavor breaking effects have already been observed. For example, experiment finds pion and kaon charge radii of r π + = 0.672 ± 0.008 fm and r K + = 0.560 ± 0.031 fm [1], which reveals flavor breaking effects of around 10%. More striking perhaps is for the neutral pion and kaon, where r π 0 vanishes and r K 0 = −0.277 ± 0.018 fm [1]. Similar effects are expected in the pion and kaon parton distribution functions (PDFs). However, while some data exists for the pion, from pion induced Drell-Yan [3], knowledge of the kaon is even more limited with only some early data onū K − (x)/ū π − (x) [4]. The naive expectation based on quark mass effects is that s K (x) will have more support at large x (harder) while u K (x) will be concentrated at smaller x (softer). Similarly, flavor breaking would imply s K (x) is harder and u K (x) is softer than the u/d quark distributions in the pion. Existing data on the pion and kaon are not in contradiction to these naive expectations, however, the error bars are large so definitive conclusions cannot be made.
Calculations of the pion and kaon are also limited within lattice QCD as compared to the proton. Among the first calculations for the pion are for the moments in the quenched approximation [5], which were later improved [6,7]. It is only recently that the first calculation of moments using local operators was performed [8]. Lattice calculations of the x-dependence of the pion and kaon PDFs became available in the last few years [9][10][11][12][13][14][15][16] using methods like the quasi-PDFs [17,18], pseudo-Ioffe-time-distributions (ITD) [19], and current-current correlators [20][21][22]. For a recent review on these approaches see Refs. [23,24]. Using these methods, one can integrate and extract the n th moments.
In this work, we calculate the non-trival moments of the pion and kaon quark PDFs up to x 3 using lattice QCD, and explore the size of quark mass effects by comparing moments in the pion and kaon, and between the light and strange quarks in the kaon. These moments are determined by directly evaluating the associated operators, using one ensemble of gauge configurations with two degenerate light quarks, and strange and charm quarks (N f = 2 + 1 + 1) of maximally twisted mass fermions with a clover improvement. We avoid operator mixing in the x 3 moments by using three different spatial and a temporal component for the associated operator. The computation of these moments provide important insight into the large-x behavior of the pion and kaon PDFs, and provide a check of systematic errors associated with recent methods to determine the full x-dependence of the PDFs using, e.g., quasi-PDFs [15,17], pseudo-PDFs [25,26], and current-current correlation methods [16,22]. With x 3 at hand for the pion and kaon, and the lower two moments calculated in Ref. [8], we make an attempt to reconstruct the x-dependence of the PDFs by fitting to a functional form, which is constrained by the baryon number sum rule. Sensitivity to the exponent β of the reconstructed PDFs q(x) (1 − x) β near x ∼ 1 is explored. Furthermore, we discuss the size of SU(3) flavor breaking effects by comparing the pion and kaon moments as well as the reconstructed PDFs.
The paper is organized as follows: In Sec. II, we present the theoretical framework and the decomposition to obtain the x 3 quark moments for the pion and kaon. We refer the reader to Ref. [8] for the formalism associated with the x and x 2 quark moments. In Sec. III we describe how to perform the non-perturbative renormalization of the three-derivative operator associated with x 3 and in Sec. IV we provide details on the analysis methods used to obtain the pion and kaon moments. Our results for the pion and kaon quark PDF moments are presented in Sec. V, where comparisons are made to other lattice QCD results and various model calculations. In Sec. VI, we present the reconstruction of the x-dependence of the PDFs and in Sec. VII we summarize this work.

II. THEORETICAL AND LATTICE SETUP
The meson matrix elements connected to x 3 , M |O {µνρσ} |M , with |M being a meson state, contain a bilinear fermion vector operator with three covariant derivatives, that is where ψ is a quark field. Curly brackets denote symmetrization over the indices and subtraction of the trace. In general, any of the indices µ, ν, ρ, σ can be spatial (1,2,3) or temporal (4) and in any combination. However, to avoid mixing with other operators we choose all four indices to be different from each other [27][28][29][30], which leads to the operator O {1234} . There are twenty-four permutations for the indices, which are all calculated and averaged over.
The relevant decomposition of the meson matrix element of the operator O {µνρσ} in Euclidean space leads to three generalized form factors, via the expression in a general frame with initial momentum p and final momentum p . P is defined as the average of the initial and final momenta of the meson, P = (p + p )/2, and ∆ is their difference, ∆ = p − p. The generalized form factors A ij are only dependent on the momentum transferred squared, Q 2 . C is a kinematic factor, which depends on the normalization of the meson state. In this work, we obtain C = for a general frame, where m M is the mass of meson M and E(p)= m 2 M + p 2 is the energy at momentum p. To extract the moment x 3 ≡ A 40 (0), we study the forward-kinematics limit of Eq. (2), that is, p = p = 0, Q 2 = 0, which leads to The kinematic factor C is simplified to 1/(2E). Similarly to x 2 , the kinematic coefficient of x 3 becomes zero in the rest frame ( p = p = 0) at the forward limit, unless all the indices of the operator are temporal. However, the renormalization of O 4444 is very complicated as it is not multiplicative and requires one to disentangle x 3 from with lower-dimension operators [27][28][29]. The only possibility to avoid any mixing in x 3 is to use the operator O {1234} , which we use in this work, and obtain: We employ the simplest boosted frame setup to get non-vanishing matrix elements, that is, all spatial components of the four-momentum are equal, p = p = (iE, p 1 , p 2 , p 3 ). To have the highest possible signal-to-noise ratio in this setup we choose the smallest nonzero value for the spatial components, that is p = p = 2π L (±1, ±1, ±1), where L is the spatial extend of the lattice used. To increase statistics, we use all eight permutations of the momentum boost, and all twenty-four combinations of indices that enter the symmetrization of O {1234} .
We calculate the connected contributions to x 3 , shown in the pictorial representation of Fig. 1. We use an ensemble of gauge configurations labeled cA211.30.32, which has been produced by the Extended Twisted Mass Collaboration (ETMC) [31]. The ensemble uses N f = 2 + 1 + 1 twisted clover fermions with clover improvement and the Iwasaki improved gluon action. The fermionic action is written as the sum of the degenerate light (S tm ) and non-degenerate heavy (S h tm ) fermion actions with: Here, χ = (u, d) represents the light quark doublet and χ h = (s, c) the heavy quark doublet. µ is the twisted quark mass of the degenerate light doublet and µ δ and µ σ are the twisted quark masses of the heavy doublet. m and m h are the (untwisted) Wilson quark masses and D W is the massless Wilson-Dirac operator. The bare untwisted Wilson masses are tuned to the critical value m = m h = m crit , which gives automatic O(a) improvement [32] and requires no further operator level improvements. The clover term multiplied by the Sheikoleslami-Wohlert improvement coefficient c sw reduces isospin symmetry breaking effects [33], since O(a) improvement is already achieved from tuning to the critical Wilson quark mass. The key parameters of the ensemble are collected in Table I and the remaining parameters are κ crit = 1/(2am crit + 8) = 0.1400645, c sw = 1.74, aµ = 0.003, aµ σ = 0.1408, and aµ δ = 0.1521.
As in Ref. [8] in which the lower two moments of the Kaon PDF were computed, in the valence sector we use the so-called Osterwalder Seiler fermions for the strange quarks, with the same value of the bare strange quark mass. The matrix elements M (p)|O {µνρσ} |M (p) use the interpolating fields of π + and K + , that is, J π + = dγ 5 u and J K + = sγ 5 u. For the pion, we only need to calculate the up quark contribution to x 3 , as the equality G u (x, x ) = γ 5 G † d (x , x)γ 5 holds for twisted mass fermion propagators. The interpolating fields at the source and the sink are smeared using Gaussian smearing. Further details can be found in Ref. [8]. The three-point correlation functions in the forward limit are given by where t i , t, t s denote the source, insertion and sink Euclidean times, respectively. Similarly, the spatial coordinates of the source, current insertion and sink are x i , x, x s . We set the source to be at t i = 0, so that the source-sink separation is t s − t i ≡ t s . In the results presented in this paper, we focus on the u + contribution to the pion, where q + ≡ q +q. To get the total connected contribution one may use the relation x 3 u + +d + π = 2 x 3 u + π . The calculation requires a large number of statistics to control the gauge noise introduced by the covariant derivatives. The signal-to-noise ratio is also suppressed due to the use of a boosted frame, which is necessary to obtain x 3 directly in the forward kinematic limit while avoiding operator mixing. We use the smallest possible values, that is, momenta of the class p i = 2π L (±1, ±1, ±1) with p i 2 = 12π 2 L 2 . This leads to a factor of eight more computational cost, which is compensated by the reduction of statistical uncertainties by a factor of about 1/ √ 8.
We analyze 122 configurations, separated by 20 trajectories to reduce auto-correlation effects. In order to control the gauge noise, we calculate the matrix elements for x 3 for more than double the statistics as compared to our previous work on x and x 2 [8]. Based on the conclusions of Ref. [8], we use t s /a = 14, 16, 18 to reliably suppress excited states. The pion, being a lighter meson, suffers from higher statistical noise. Therefore, we add t s /a = 12 in the procedure for eliminating excited-states contamination. The statistics used for each value of t s are listed in Table II.

III. RENORMALIZATION
The three-derivative operator, in its general case, exhibits mixing with lower dimension operators [27][28][29][30]. As mentioned in the previous section, we choose all indices of the operator to be different (µ = ν = ρ = σ = µ), which avoids such a mixing and, therefore, its renormalization is multiplicative. We calculate the corresponding renormalization function, Z vDDD , non-perturbatively following the procedure we developed in Refs. [34][35][36]. We start by writing the bare vertex function as with J = γ µ D ν D ρ D σ , and u and d representing quark fields in the physical basis. p is the vertex momentum allowed by the boundary conditions, and V is the lattice volume. The Dirac and color indices of G(p) are suppressed for simplicity. We employ the momentum source approach, introduced in Ref. [37], which uses directly Eq. (8) with a source that is momentum dependent. For twisted mass fermions, we make use of the symmetry S u (x, y) = γ 5 S d † (y, x)γ 5 between the u− and d−quark propagators, and therefore, extract the vertex function with a single inversion per momentum. While this method requires separate inversions for each momentum employed in the calculation, it has the advantage of high statistical accuracy and the evaluation of the vertex for any operator at no significant additional computational cost.
Z vDDD is calculated by employing the Rome-Southampton method (RI scheme) [38], which involves the amputated vertex function S u (p) and S d (p) are the propagators in momentum space in the physical basis, defined by In practice, we work in the twisted basis at maximal twist, in which Eq. (8) takes the form x,y,z,z The above expression can be simplified by using S u (x, z) = γ 5 S d † (z, x)γ 5 . The renormalization function in the RI scheme are determined by the conditions where the trace is taken over spin and color. The momentum of the vertex function is indicated by p, and is set to the RI renormalization scale, µ 0 . S Born (Γ Born vDDD ) is the tree-level value of the fermion propagator (operator). We differentiate Z vDDD and Z vDDD , as the former depends on the pion mass of the ensemble and the initial RI renormalization scale. The latter is our final estimate after the chiral extrapolation and after the limit (aµ 0 ) 2 → 0 has been taken.
We evaluate the vertex functions and propagators for a wide range of values for (a p) 2 , using 10 gauge configurations, which leads to per mil statistical accuracy. We use momenta that have the same spatial components, that is: [2,9] , n x [2,5] , (ap) 2 ∈ [0.9 − 6.7] , where L t (L s ) is the temporal (spatial) extent of the lattice. Democratic momenta in the spatial directions that satisfy, in addition, Eq. (14) reduce non-Lorentz invariant contributions . This is based on empirical arguments [39] and is being implemented in all calculations by our group. We further improve Z vDDD by subtracting the O(g 2 a ∞ ) artifacts from Z q , which enters the renormalization condition of Eq. (12). The artifacts are computed to one loop in perturbation theory and to all orders in the lattice spacing, O(g 2 a ∞ ), as outlined in Refs. [36,40]. Note that the vertex function of the three-derivative operator also contain O(g 2 a ∞ ), but have not been calculated yet.
We obtain Z vDDD on five N f = 4 ensembles at the same lattice spacing as the N f = 2 + 1 + 1 ensemble we use for the meson matrix elements. The ensembles correspond to different pion masses, and are used in order to take the chiral limit. The parameters of the N f = 4 ensembles used for Z vDDD are given in Table III. The chiral limit is taken using a quadratic fit with respect to the pion mass, or linear in aµ, giving a zero slope within uncertainties in both cases. The chirally extrapolated values for Z vDDD in the RI scheme are converted and evolved to MS(2 GeV) using an intermediate Renormalization Group Invariant (RGI) scheme. Finally, a linear fit with respect to (aµ 0 ) 2 is applied to the MS estimates to eliminate residual dependence on µ 0 , that is Z O corresponds to the final value of the renormalization function for operator O. The estimates for Z vDDD in the RI and MS schemes as a function of the initial RI renormalization scale are shown in Fig. 2. Z MS vDDD is given at a scale of µ = 2 GeV. We find that, for (aµ 0 ) 2 ≥ 3 the purely non-perturbative data exhibit a small residual dependence on the initial scale µ 0 they were evolved from. The subtraction of the lattice artifacts in Z q results in a smaller slope, demonstrating the effectiveness of the artifact-subtraction procedure. We obtain with the number in the first (second) parenthesis being the statistical (systematic) uncertainty. The indicated systematic effect is obtained by taking the deference between the estimates for the fit interval (a µ 0

A. Setup
To extract the meson matrix elements, one requires to take the ratio of the three-point correlation functions of Eq. (7) with the two-point functions, given by The normalization of the meson state is 0|J The ground state contribution can be isolated from the ratio which cancels unknown overlap terms between the interpolating field and the meson state. In this work, we apply two methods to obtain the ground state, namely a single-state (plateau) fit, and a two-state fit.

Two-point function
In both methods for identifying and eliminating excited states, analyse the two-point function and perform a single-state and two-state fit, the latter being In the above equation, the fitting parameters are the ground state energy, E 0 (p 2 ), the first excited state energy, E 1 (p 2 ), and the amplitudes c 0 and c 1 . The plateau fit takes into account only the first term of Eq. (19). We fit the two-point functions averaged over the eight different directions of meson momentum (p 2 = 12π 2 /L 2 ), as the energies only depend on the momentum squared. In addition, using the averaged two-point functions improves the stability of the fit compared to fits on the individual momentum directions. The effective energy is calculated using the formula which assumes that the two-point functions are symmetrized, i.e., that the two-point functions at t have been averaged with the values at T − t. For the fits, we select a range of t ∈ [t low , 31] with varying t low . We choose the lowest values of t low , such that the resulting ground state energies, E 2−state and E plat satisfy the condition where δE plat is the error on the plateau fit. This criterion, while not unique, works very well on the lattice data we obtain here (see Fig. 3 and Fig. 5).

Single-state fit
We calculate the ratios in Eq. (18) using a fit of the two-point function instead of the actual lattice data, that is where c 0 and E 0 are the ground state amplitude and energy, respectively. Using the modified two-point function removes t s dependence from the plateau values of the the ratios so that the plateaus converge at high t s . Based on the plateau method, at insertion times far enough from the source and sink, the above ratio becomes time independent, i.e., In practice, we apply a constant fit in a region where a plateau is identified. The time-independent ratio (plateau) is renormalized multiplicative with Z vDDD , and is related to the desired x 3 as given in Eq. (4) Indeed, we confirm that the signal is found in the imaginary part of the three-point function. The single-state fit is applied on each t s separately, and the ground state is extracted at the lowest t s , beyond which the plateau value is t s -independent. We also compare with the two-state fit to confirm convergence.

Two-state fit
For the two-state fit we include the data for all t s simultaneously. The two-state fit is calculated by fitting the three-point functions to the Ansatz where the fitting parameters are the amplitudes A 00 , A 10 , and A 11 since, at zero momentum transfer, A 10 = A 01 .
To avoid heavy notation we do not include a subscript M in the parameters and energies. We use the energies E 0 and E 1 calculated from the two-state fit on the two-point functions. The results of the two-state fits on the two-and three-point functions are related to the matrix elements by so that x 3 M is calculated from the two-state fit as As mentioned above, the two-state fit is useful to confirm ground-state dominance beyond a certain value of t s .

B. Pion
We start our presentation with the extraction of the ground state energy for a boosted pion, which is needed for the analysis of the three-point functions. We follow the setup outlined above, and the results are shown in Fig. 3. We apply a single-state (plateau), as well as two-state fits, and test the results against energies obtained via the continuum dispersion relation E 2 (p) = m 2 + ( 2π L p) 2 , where m is obtained from the effective mass of the two-point correlation function at zero momentum, i.e. Eq. (20). For the two-state fit we vary t low between t = 1a and t = 7a. The single-state fit is applied for t low ∈ [5a − 14a]. As can be seen in Fig. 3, we find that the ground-state energy is isolated already at t low /a = 2 for the two-state fit and at t low /a = 8 for the plateau fit. The ground-state energy obtained from the plateau fit is aE plat = 0.3668 (24) and from the two-state fit is aE 2−state = 0.3674 (18). These values satisfy the criterion of Eq. (21), and also, are compatible with each other and with the dispersion relation.
Using the above values for the ground state energy we form the modified two-point function of Eq. (19) and take the ratio of Eq. (18). We also ensure that excited-state effects can be successfully suppressed by calculating the threepoint correlators at different source-sink time separations. In our previous study [8], we calculated x in the rest frame for t s /a = 12, 14, 16, 18, 20, 24 and found that the results converge at t s /a ≥ 18. Therefore, for the boosted frame we focus on t s /a = 12, 14, 16, 18. Fig. 4 shows the ratios which lead to x 3 u π , that is, the right-hand-side (rhs) of Eq. (24) for the plateau method and the rhs of Eq. (27) for the two-state fit. We find that the signal is more noisy for t s /a ≥ 16, due to the three covariant derivatives and the boosted frame setup. All plateau fits, however, are found to be consistent with one another as well as with the two-state fit result. This finding suggests that any remaining excited-state effects are within the reported errors. The plateau and two-state fit values are collected in Table IV. For completeness, we give the updated values for x 2 u π obtained with double the statistics for t s /a = 14, 16, 18 compared to our previous work [8]. The computational challenges associated with the gauge noise contamination in three-derivative operators, the boosted frame with three nonzero spatial components, coupled with the light mass of the pion are reflected in the increased uncertainties in x 3 u π . However, this is not the case for the kaon, which is about twice heavier than the pion for this ensemble (see Table V).

C. Kaon
Our analysis for the kaon follows the same procedure for the two-point and three-point correlation functions, as in the case of the pion. In Fig. 5, we show the extraction of the ground-state energy, E K . The dependence of E K on t low is similar to the pion. We extract the plateau value aE plat = 0.4230(12) from t low = 11a, and the two-state fit value aE 2−state = 0.4230 (7) at t low = 1a. In Fig. 6, we plot the ratio leading to x 3 K for the up and strange quark. We show the four values of t s and compare with the two-state fit. As can be seen, the gauge noise is decreased compared to the pion, attributed to the heavier mass of the kaon. There is a clear signal for both flavors of the kaon for all source-sink time separations. Similarly to the pion, the plateau values are consistent with the two-state fit. The plateau and two-state fits are collected in Table V.  The ratios leading to x 3 for the kaon. The top and bottom panels correspond to the up and strange contributions, respectively. The notation is the same as in Fig. 4.

A. Results
Based on our analysis, we choose the results extracted using two-state fits as our final values for all quantities. We report as systematic uncertainties the difference between the value extracted using the two-state fit and that determined using the plateau fit at t s = 18a. Our final results in the MS at a scale of 2 GeV are for the second moment, and for the third. We also calculate the ratios x 3 / x and x 3 / x 2 . The calculation of x and x 2 is discussed in Ref. [8]. These ratios are of interest because they reveal partial information on the x-dependence PDFs, in particular, the large-x region. In general, the statistical errors for x 3 are larger than the lower moments, which is propagated to the ratios. The x 3 to x ratios are and the x 3 to x 2 ratios are It is also interesting to compare x n π and x n K , as it is related to the SU(3) flavor symmetry breaking due to the heavier mass of the strange quark. This is important because some quantities may be sensitive to SU(3) flavor symmetry breaking, such as the pion and kaon radii [1]. For the comparison of the up quark contributions we find while for the strange quark in the kaon over the up quark in the pion we have We find that the SU(3) symmetry breaking is ∼ 5 − 10% for x and ∼ 10 − 20% for x 2 . The results for x 3 indicate a symmetry breaking of ∼ 30 − 50%, with larger uncertainties. These results are very interesting because based on intuitive arguments, the strange quark PDF has its support at higher x than the up quark PDF, which indicates that the symmetry breaking is more pronounced in the higher moments. We will discuss further the SU(3) flavor symmetry breaking in Section VI.

B. Other lattice calculations
There are only a few direct calculations for x 3 extracted from the three-derivative local vector operator, starting with the pioneering work of QCDSF-UKQCD in 1997 [5], which was later extended in 2007 [6,7], and reanalyzed in Ref. [41]. All these calculations are for the pion; our results are the first reported for the kaon, using the threederivative local operator. Interest in the pion and kaon structure has recently been renewed, due to novel approaches to extract the x-dependence of PDFs, such as the quasi-PDFs [17,18], the pseudo-Ioffe-time-distributions (ITD) [19], and current-current correlators [20][21][22]. Using such methods, the lowest moments of the pion [12,14,15], and the kaon [14] have been obtained [12,14,15], either from integration on the x-dependent PDF, or via the so-called operator product expansion without an operator product expansion (OPE without OPE) method [42]. Given the small number of calculations, we compare our results with all these methods. One has to bear in mind that, each calculation has its own systematic uncertainties and uses a different methodology. Therefore, the comparison is qualitative at this stage, as not all sources of systematic uncertainties have been quantified.
The calculation of Ref. [5] was done in the quenched approximation for Wilson fermions using ensembles with pion mass 712, 1013, 1208 MeV. They used an operator with only two indices different (O v4 = O ii44 ), which requires momentum boost only in one spatial direction. However, this operator exhibits mixing with lower dimensional operators, which is difficult to eliminate. The calculation ignored the mixing and x 3 π was renormalized using results from perturbation theory. The reported values for x 3 u π in the MS scheme at a scale of 2.4 GeV are 0.0619(45) 0.0580(65) 0.054 (18) for the ensembles with pion mass 1208, 1013, 712 MeV, respectively. In Table VI we present their estimate for the extrapolated value, 0.048 (20), after evolution to 2 GeV. The results of Ref. [5] have been analyzed in Ref. [41] using different methods to perform the chiral extrapolation on the three ensembles mentioned above. For two of the methods they report x 3 u π = 0.043 (15)(3), and 0.05(2) at a scale of 2.4 GeV. The value of the first method is evolved to 2 GeV and is given in Table VI.
Almost a decade after their first calculation, QCDSF-UKQCD has improved their work in more than one ways, as presented in Refs. [6]. The ensembles employed are unquenched (N f = 2) clover fermions. Several ensembles were used to extract x 3 π with pion mass between 450 -1180 MeV. The operator O v4 was used, and while the mixing was ignored, the renormalization was done non-perturbatively. Their preliminary results chirally extrapolated to the physical pion mass gave x 3 π = 0.074(9)(4) in the MS at a scale of 2 GeV. It should be noted, that the accuracy of the chirally extrapolated value is heavily influenced from the accuracy of the ensembles with m π = 800 MeV and higher. Another update of the calculation was presented in the thesis of Ref. [7], using N f = 2 clover fermions at four β values, several lattice spacings between 0.068 -0.115 fm and a wide range of pion mass values between 440 -1173 MeV. They report x 3 u π =0.074(9)(4) in the MS scheme at 2 GeV at the physical pion mass obtained from a chiral extrapolation.
Much more recently, there have been explorations of the x-dependent pion and kaon PDFs [9][10][11][12][13][14][15][16], and some extract the first moments with indirect methods. Unlike our work, and the aforementioned calculations, these new methods calculate matrix elements of non-local operators, with the quark fields separated by a finite spatial distance connected through a straight Wilson line. Therefore, the mixing observed in the moments of PDFs from towers of n-derivative local operators, is not relevant here. However, these methods have other systematic uncertainties, due to the need of boosted meson states, namely either requiring the boost to be large or the product of spatial separation and boost to be large. Furthermore, a matching kernel is necessary to relate these matrix elements to the light-cone PDFs.
in Ref. [12] two ensembles of N f = 2 + 1 clover fermions with pion mass 415 MeV and two volumes were used. They extracted the x-dependence of PDFs using the pseudo-ITD approach. They used two methods to extract the lowest moments: OPE without an OPE and an integration of the PDFs that gives x 3 u π = 0.046 (19) in the MS at 2 GeV. Ref. [14] employs the quasi-PDFs method using a mixed-action setup (clover on HISQ). The calculation was performed on three ensembles with pion mass 217, 310, 319 MeV and two lattice spacings (a = 0.06, 0.12 fm). x 3 u π was renormalized non-perturbatively and a value of x 3 u π = 0.057(10) at a scale of 5.2 GeV is reported after chiral extrapolation to the physical point. A similar analysis for the kaon gives x 3 u K = 0.042(6) and x 3 u K = 0.070 (6), but the scale is not reported. The result for the pion is presented in Table VI after evolution to 2 GeV. Finally, the work of Ref. [15] explored both the quasi-PDFs and pseudo-ITD approaches, and also used a mixed-action setup of clover valence fermions on N f = 2 + 1 HISQ configurations. The pion PDF was calculated using two ensembles at m π = 300 MeV and lattice spacing 0.04 and 0.06 fm. The third moment is obtained with fits to the pion PDF with a 2-and a 4-parameter Ansatz. They report the valence case x 3 π = 0.0652(49)(36), 0.0647(47) (38), for the 2-parameter and 4-parameter fits, respectively. The results are given in the MS scheme at 3.2 GeV. We extract x 3 π at 2 GeV using the finest lattice and the 2-parameter fit, which is provided in Table VI. For completeness, we summarize our calculation for the kaon in  [7] local operator non-perturb. present chiral extrap. 2 0.074(10) 2 GeV Ref. [12] pseudo-ITD non-perturb. N/A 415 2+1 0.046 (19) 2 GeV Ref. [14] quasi-PDF non-perturb. N/A chiral extrap. 2+1+1 0.073(13) 5.2 GeV Ref. [15] pseudo-ITD non-perturb. N/A 300 2+1 0.075 (61) 3.2 GeV TABLE VI. Comparison of lattice results for x 3 u π in the MS scheme at 2 GeV. The evolution from the reported scale ("initial" scale) to 2 GeV is applied to NNLO. Statistical and systematic uncertainties have been added in quadrature where applicable.
As can be seen in Table VI, there is a range of values obtained with different methods. We find that, our results are compatible within uncertainties with the calculations of Refs. [5,41], which, however, used an operator that exhibits mixing. Ref. [5] employs a perturbative renormalization prescription, which justifies the higher value compared to the other studies. The value of Refs. [7] is higher than the other calculations with local operators. The comparison with the indirect methods to extract the Mellin moments from integration or fits on the pion PDF, also shows compatibility with our value, except the one for Ref. [14], which is at the high end. It should be noted that these calculations carry very large uncertainties (∼ 15% − 95%) and the comparison is inconclusive. From these results one can extract the range in which x 3 π is. It is worth mentioning that our value is at the lower end of the range, which is a consequence of the suppression of excited states. For example, our results for t s /a = 14, which is at the range used in other calculations is higher (0.031 (15)). We emphasize that we obtain directly x 3 using a local operator that avoids mixing with lower-dimensional operators, but requires all spatial components of the meson momentum boost to be nonzero. The lattice data are renormalized non-perturbatively, and the renormalization function is multiplicative.

C. Model calculations and global fits
There are a few model calculations and global fits on experimental data for the three-four lowest moments of pion and kaon PDFs, and we compare these here with our results. We emphasize that the comparison is only qualitative, as many of the calculations do not have quantified uncertainties. Also, our calculation is at higher-than-physical pion and kaon masses and only the connected diagram is included.
One of the first calculations is a next-to-leading-order analysis of several π ± N experimental data, including Drell-Yan and prompt photon production, presented in Ref. [43]. They obtain x 3 u π = 0.058(4) at a scale of 2 GeV. Much later, an updated analysis of the moments of pion PDFs to next-to-leading-order using the Fermilab E-615 pionic Drell-Yan data was carried out and can be found in Ref. [44]. A value of x 3 u π = 0.045(3) is given at 5.2 GeV. The JAM global fit analysis is performed at 1.3 GeV for the third non-trivial moment of the pion [45], and the obtained value is x 3 u π = 0.074. Ref. [46] presents a calculation of the valence quark PDF for the pion using Schwinger-Dyson equations (DSE) and obtains x 3 u π = 0.049 at 2 GeV. A more recent DSE study can be found in Ref. [47] with x 3 u π = 0.052 at 2 GeV. Ref. [48] reports x 3 u π = 0.049(7) at 2 GeV using the Bethe-Salpeter equation (BSE). The recent calculation of Ref. [49] applied a rainbow-ladder truncation of DSEs, and therefore all planar diagrams were summed and the non-perturbative gluon dressing of the quarks was correctly accounted for. They find x 3 u π = 0.109 at 0.78 GeV. A calculation using the chiral constituent quark model is found in Ref. [50], which gives x 3 u π = 0.048 at a scale of 5.2 GeV. Finally, Ref. [51] combined QCD evolution with light front quantization to obtain the pion PDFs, as well as the moments up to x 4 . At 2 GeV they find x 3 u π = 0.057 (8). The kaon x 3 was also studied in Ref. [47] and the results are x 3 u K = 0.048 and x 3 s K = 0.092 at 2 GeV. The model calculation of Ref. [50] gives x 3 u K = 0.045 and x 3 s K = 0.049 at 5.2 GeV. The findings of Ref. [49] are x 3 u K = 0.092 and x 3 s K = 0.143 at 0.78 GeV. Ref. [51] report x 3 u K = 0.050(6) and x 3 s K = 0.066(9) at 2 GeV. The aforementioned results can be found in Table VIII for the pion and Table IX for the kaon, after evolution to 2 GeV. An extended list can be found in Ref. [51]. We find that our data for the pion and up part of kaon are lower than most of the other calculations. However, some of the calculations do not include systematic uncertainties, which prevents a meaningful comparison. A better agreement is observed for the strange part of the kaon with our values being in the middle of the range from the other calculations. Let us remind the reader that our calculation, like all other lattice results mentioned above, focuses in the connected diagram. Also, the ensemble used has a pion mass of 260 MeV and a kaon of 530 MeV, which are higher than the physical values.
Reference and x 3 s K in the MS scheme at 2 GeV. The evolution from the reported scale ("initial" scale) to 2 GeV is applied to NNLO. Statistical and systematic uncertainties have been added in quadrature where applicable.

VI. RECONSTRUCTION OF PDFS
It is generally believed that the reconstruction of PDFs from their Mellin moments is, at best, challenging on the lattice for a number of reasons. The signal-to-noise ratio decays fast with the addition of derivatives in the operator, requiring increased computational cost so that gauge noise is controlled. The moments x 2 and x 3 can only be obtained in a kinematic framework where the hadron has momentum with nonzero spatial components. In fact, to avoid mixing with lower dimensional operators under renormalization, the initial and final states should carry momentum with at least two and three nonzero spatial components, respectively. Such a setup comes at increased computational cost. In addition, the mixing under renormalization for the moments with n > 3 cannot be avoided, regardless of the kinematic framework. Because of these challenges, early attempts have been inconclusive in determining whether it is feasible to reconstruct PDFs from lower moments using lattice QCD (see, e.g., Refs. [52,53]). More recently, methods to extract higher moments have been proposed, using smeared operators [54], heavy-quark operator product expansion (HOPE) [55,56], and light-quark current-current correlators [57].
While the above challenges are true, numerical simulations have advanced significantly with more computational power, better algorithms, implementation of non-perturbative renormalization, and methods to control gauge noise. In fact, in this work and in Ref. [8], we demonstrate that the Mellin moments with n < 4 can be obtained with reliable elimination of excited states and in a setup that does not contain mixing with lower dimensional operators. Therefore, we attempt the reconstruction of the x-dependence of PDFs. The goal is threefold: a. understand the limitations of the reconstruction; b. study the large-x behavior; and c. extract the moments with n > 3 from the reconstructed PDFs.

A. Setup
We use the standard functional form to obtain the x dependence of the pion and kaon PDFs, q f M (x), where (M, f ) = (π, u), (K, u), (K, s). N is a normalization defined by charge conservation leading to where B is the Euler beta-function. The fit parameters in Eqs. (46), (48) are α, β, γ and ρ. Their values depend on M and f , but we omit the subscript and superscript in the following equations for simplicity in the notation. The parameter ρ is generally assumed to be negligible [47], and therefore, we omit the term ρ √ x. By integrating Eq. (46), we can extract the n th -moment as a function of the fit parameters, that is Our results for x n , n = 1, 2, 3 are used as input for Eq. (49) to extract the fit parameters. The results of x are given in Eqs. (38) -(40) of Ref. [8]. For x 2 we use the values given in Eqs. (28) - (30), which have been obtained at higher statistics than the values reported in Ref. [8]. x 3 is given in Eqs. (31) - (33). To compare with results from global fits and models, we evolve our results for the moments to a scale of 5.2 GeV using NNLO expressions.

B. Lattice data
We apply a 2-and a 3-parameter fit to examine the effects on the PDF reconstruction. In the case of the 2parameter fit we set γ = 0, as used in many such fits. The values extracted for the parameters from the fits are given in Table X for the pion and kaon. We find that the fits for the pion are less stable than the ones for the kaon, due to the enhanced gauge noise in the former. For the case of the pion and the strange contribution of the kaon, we find that the parameter γ has large uncertainties. This is due to the fact that the number of moments are not enough to perform a 3-parameter fit, and the parameters have a competitive role in the fit. For instance, we find that in the pion and strange kaon, γ = 0 within errors.
Using the parameters of Table X obtained from both fits we reconstruct the x-dependent PDFs as shown in Fig. 7. We find that the shape of the PDFs has a mild dependence on the choice of fit, mostly in q u K . The uncertainties fit type α u π β u π γ u π χ 2 /d.o.f.  for q s K are increased for the 3-parameter fit, while they are very similar for the other two PDFs. For the remaining presentation we focus on the 2-parameter fits.
Next, we want to study the effects of excited-states contamination on q f M (x) that may be non-trivial because the dependence of the fit parameters in Eq. (46) on the moments is non-linear. We apply the 2-parameter fit on our results at t s /a = 14 − 18, as well as the two-state fit and reconstruct the PDFs, as shown in Fig. 8. We observe a nice convergence for all the cases as t s increases. Similarly to the behavior of the moments, we find that the excited-states contamination at t s 1.4 fm leads to a PDF that is higher than the two-state fit and the plateau fit at t s 1.5 fm. Eventually, the two-state fit values converge with a peak around x ∼ 0.3 − 0.4, that is xq u π (0.3) ∼ 0.4, xq u K (0.3) ∼ 0.4 and xq s K (0.4) ∼ 0.5. We use as final PDFs those extracted from a 2-parameter fit on the two-state fit results for the moments (purple band in Fig. 8).
One concern is whether the use of the n ≤ 3 moments can successfully reconstruct the PDF. To address this question we use the JAM data for the pion PDF [45], as well as its moments. We follow the same procedure as for our lattice data, that is, fit the parameters α, β, γ using as input only the JAM moments with n ≤ 3, and then produce the reconstructed PDF via Eq. (46) with ρ = 0. The reconstructed PDF is then compared to the original JAM PDF and is shown in Fig. 9. We use bootstrap sampling to obtain the uncertainties and width of the band. As can be seen, the two are in very good agreement within uncertainties for almost all x region. The reconstruction of the PDF from its moments with only n ≤ 3 leads to much larger uncertainties due to the truncation of the information that is used to extract the PDF. Furthermore, we use the reconstructed PDF to estimate the n = 4 moment via Eq. (49), and we find x 4 u π = 0.026 (2). This is in excellent agreement with the moment as extracted from the JAM framework, x 4 u π = 0.027 (2), which is more accurate, as expected. These results suggest that the reconstruction of the PDFs using the n ≤ 3 moments is indeed feasible. To further test the sensitivity of our fits, we change the number of inputs used. In particular, we perform fits including moments up to x nmax , with n max =2, 3, or 4. For n max =2 and 3 we only use lattice data, while for n max =4 we add another constraint by using the value of x 4 from global fits and models for the pion and kaon, respectively. In particular, we use x 4 u π = 0.027(2) from the JAM analysis [45], and x 4 s K = 0.029 +0.005 −0.004 , x 4 u K = 0.021 +0.003 −0.003 from BLFQ-NJL [51]. It should be emphasized that, combining lattice data with results from model calculations is a useful exercise for understanding the effect on the fits. However, it is not a preferred direction due to the various sources of uncertainties, and only lattice data enter our final results. We remind the reader that our calculation is performed at a pion mass of 260 MeV and a kaon of 530 MeV, which are larger than their physical values. While for the kaon this is only 7% larger, for the pion this corresponds to a factor of about two. However, it was shown that the pion moment x has insignificant pion mass dependence [58], which implies the same for the higher moments. Therefore, one can neglect the pion mass dependence when combining our results with the JAM moment x 4 . In Fig. 10 we compare the resulting PDF ( x nmax = x 4 ) to the one that uses only lattice data ( x nmax = x 2 , x 3 ). The results show that the addition of n = 3 improves the constraint of PDFs. Interestingly, the addition of n = 4 does not affect the shape of the PDFs. Therefore, the effect of higher moments is within the shown uncertainties. To summarize, our final estimates are obtained from a 2-parameter fit applied on the two-state fit results using our lattice results up to The large-x behavior of the pion and kaon PDFs has been of great interest, due to different findings between existing data and model calculations. For the pion, the analysis of the pion Drell-Yan data from the Fermilab E615 experiment [3] suggests a (1 − x) 1 fall (β = 1), while the data of Ref. [59] indicate a (1 − x) 2 dependence (β = 2). Also, DSE results [47] find a coefficient β closer to 2, in support of the arguments that the distribution at large-x is dominated by the term (1 − x) 2+γ , where γ ≥ 0 is an anomalous dimension. One can argue either direction: Ref. [3] uses evolution equations to only leading order in perturbative QCD, which is not sufficient for convergence. On the other hand, the model calculations have limitations on how well they describe QCD. To date, the tension persists. Therefore, it is desirable to address this issue from lattice QCD. Recently, the large-x dependence was discussed using new methods to access the x-dependence of PDFs (see, e.g., Refs. [10-12, 14, 16, 60]). In this work, we address the large-x behavior using the reconstructed PDFs from their Mellin moments. As can be seen in Table X and in Fig. 7, our data show a preference in the functional form ∼ (1 − x) 2 for both the pion and kaon.
The pion and kaon PDFs can be compared in order to address SU(3) flavor symmetry, which is broken due to the larger mass of the strange quark as compared to the up and down quark. This leads to a mass difference between the pion and kaon, a manifestation of the SU(3) flavor symmetry breaking effect in the Nambu-Goldstone bosons. To this end, we compare the pion and kaon PDFs in Fig. 11. The distributions xq u π (x) and xq u K (x) are in full agreement for all regions of x, besides a minor tension around x = 0.5. Based on this behavior, one can argue that the up quark plays an equal role in the pion and kaon PDF, and has most of its support in the small-to intermediate-x regions. The strange quark in the kaon shows a tension with the up quark between x = 0.3 and x = 0.8, with the strange quark having more support in the large-x region, as one would expect from quark mass effects. We find that the peak of the distributions are xq u π (x = 0.3) = 0.43 (5), xq u π (x = 0.28) = 0.42 (2), and xq u π (x = 0.36) = 0.51 (2). Using the reconstructed pion and kaon PDFs we apply the appropriate integrals to extract their moments. One of the advantages of extracting the moments beyond n = 3 from the x-dependent PDF is that it avoids the operator mixing problem, which one has to deal with in the calculation of matrix elements using n th -derivative operators. The moments up to x 6 in the MS scheme at 27 GeV 2 are given in Table XI. Namely, we use Eq. (49) and the parameters obtained for the 2-parameter fit tabulated in Table X. We propagate both the statistical error and the systematic error of residual excited states contamination as given in the first and second parenthesis, respectively, on the mean values in Table XI. Note that the statistical uncertainties are well controlled for x n with n > 3, something that would not be feasible had these moments been calculated directly as matrix elements of the nucleon with the same statistics. Our value for x 4 is in agreement JAM value x 4 u π = 0.027 (2).

C. Comparison with other studies
As mentioned in the previous section, there are a few calculations of the pion and kaon PDFs, which we compare to our results in Fig. 12. In the left upper panel, we compare with the lattice results extracted from the pseudo-ITD approach [12] and the current-current correlators (LCS) method [10,16]. We find best agreement with the pion PDF obtained by three parameters fit using the LCS method, while the pseudo-ITD data have a lower peak. We include in the plot the E615 [3] and ASV data [59]. The latter include soft-gluon resummation and use the next-to-leading order formalism. As can be seen, while the E615 and ASV data are in agreement in the intermediate x-rang the ASV data fall off faster. The global fits of the JAM Collaboration [45,61] is also shown. The JAM fit describes well the E615 data. In the upper right panel we compare with Dyson-Schwinger (DSE) [47] as well as the updated DSE '18 [49], the chiral constituent quark model (χCQ) [50], and the BLFQ Collaboration results from in the light front quantization and QCD evolution (NJL) [51]. We note that the pion PDF can also be extracted from the determination of the nucleon PDF from light-front holographic QCD [62]. For the kaon, there exist limited calculations. In the lower panel of Fig. 12 we compare with the χCQ results [50], the BLFQ-NJL data [51], and DSE '18 [49] for the up quark (left) and strange quark (right). All results are in agreement in the small-x region (x < 0.1) for the pion and xq u K (x). The χCQ have a different slope in the small-x region. For x > 0.6 our results agree with all other results except the original E615 data. The intermediate region reveals disagreement between the various methods; our pion results agree with DSE [47] for all regions of x, but overestimate the peak as compared to DSE'18. For the kaon we find that there is qualitative agreement in the small-and large-x regions, while there is a tension in the intermediate-x region. in all cases, our dta are larger in the intermediate x region. However, we need to stress that there are no experimental data, and therefore, the comparison is qualitative, as all calculations carry non-quantified systematic uncertainties. Our results (blue band) use data up to x 3 obtained with the two-state fits analysis and a 2parameter fit. The E615 data [3] are shown with gray points, the rescaled ASV curve [59] with cyan color. The JAM global fit is shown with pink band, and the lattice results from pseudo-ITD [12] and current current correlators (LCS) [16] are shown with orange and green band, respectively. Top right panel: Comparison of our results for xq u π (x) with DSE [47] and the updated DSE'18 [49], BLFQ-NJL [51] and χCQ [50]. Bottom panel: Same as top panel for the kaon xq u K (x) (left) and xq s K (x) (right).

VII. SUMMARY
In this paper we present a calculation of the Mellin moment x 3 for the pion and kaon. We use one N f = 2 + 1 + 1 ensemble corresponding to a pion mass of 260 MeV and a kaon mass of 530 MeV. A momentum-boosted kinematical frame is required to access x 3 from Eq. (4). In fact, the momentum boost should have all spatial components nonzero to avoid mixing under renormalization. Here we use the minimum momentum, that is | p i | = √ 12π L (0.72 GeV). We find that this momentum is small enough and momentum smearing [63] does not have an advantage. Excited-states contamination are studied using four values of t s ranging from 1.12 fm to 1.67 fm. We perform a single-state and a two-state fit to ensure ground state dominance. Our analysis shows that t s > 1.5 fm is sufficient to suppress excited states.
Within this work, we also calculate non-perturbatively the renormalization function of the three-derivative operator with all Dirac indices unequal, a choice that avoids mixing. We use the RI scheme, and then convert to the MS scheme and apply evolution to 2 GeV using the RGI intermediate scheme. The results for x 3 are given in Eqs. (31) - (33) in the MS at 2 GeV. We also study ratios of moments including the moments for the kaon over pion, which are connected to SU(3) flavor symmetry breaking. We find that for the low moments this breaking is 5-10%, and increases up to about 30-50% for x 3 . Such an effect indicates that the strange quark in the kaon has its support at higher x values than the up quark in the pion and kaon.
One of the interesting aspects of this work is the reconstruction of the pion and kaon PDFs from the Mellin moments up to x 3 , as described in Section III. We apply 2-parameter and 3-parameter fits using the standard functional form for PDFs. We examine excited-states contamination and find that a two-state fit is necessary to suppress excited state effects which are not negligible for source-sink time separation below 1.5 fm. Furthermore, we reconstruct the PDFs varying the highest moment used as input to be x 2 , x 3 , or x 4 . We find that, including x 4 does not improve the PDF reconstruction.
In conclusion, we find that a 2-parameter fit of the form x α (1 − x) β using the moments up to x 3 is sufficient to reconstruct the PDF. We utilize our results on the PDFs in more than one ways. First, there has been a great interest on the value of β which captures the large-x behavior. Our lattice data exhibit a behavior of ∼ (1 − x) 2 . Second, having the functional form of the PDFs, we can obtain the higher moments, and provide results up to x 6 . Last, but very importantly, we examine the SU(3) flavor symmetry breaking in the PDFs, and find that the conclusions are consistent to those drawn when comparing ratios of the Mellin moments as described above. In particular, we see that the up quark has approximately the same contribution in the pion and kaon and that the strange quark has a more prominent contribution in the intermediate-to large-x region.