The Mellin moments $\langle x \rangle$ and $\langle x^2 \rangle$ for the pion and kaon from lattice QCD

We present a calculation of the pion quark momentum fraction, $\langle x \rangle$, and its third Mellin moment $\langle x^2 \rangle$. We also obtain, for the first time, $\langle x \rangle$ and $\langle x^2 \rangle$ for the kaon. We use an ensemble of two degenerate light, a strange and a charm quark ($N_f=2+1+1$) of maximally twisted mass fermions with clover improvement. The quark masses are chosen so that they reproduce a pion mass of about $260$ MeV, and a kaon mass of 530 MeV. The lattice spacing of the ensemble is 0.093 fm and the lattice has a spatial extent of 3~fm. We analyze several values of the source-sink time separation within the range of $1.12-2.23$~fm to study and eliminate excited-states contributions. The necessary renormalization functions are calculated non-perturbatively in the RI$'$ scheme, and are converted to the $\overline{\rm MS}$ scheme at a scale of 2 GeV. The final values for the momentum fraction are $\langle x \rangle^\pi_{u^+}=0.261(3)_{\rm stat}(6)_{\rm syst}$, $\langle x \rangle^K_{u^+}=0.246(2)_{\rm stat}(2)_{\rm syst}$, and $\langle x \rangle^K_{s^+}=0.317(2)_{\rm stat}(1)_{\rm syst}$. For the third Mellin moments we find $\langle x^2 \rangle^\pi_{u^+}=0.082(21)_{\rm stat}(17)_{\rm syst}$, $\langle x^2 \rangle^K_{u^+}=0.093(5)_{\rm stat}(3)_{\rm syst}$, and $\langle x^2 \rangle^K_{s^+}=0.134(5)_{\rm stat}(2)_{\rm syst}$. The reported systematic uncertainties are due to excited-state contamination. We also give the ratio $\langle x^2 \rangle/\langle x \rangle$ which is an indication of how quickly the PDFs lose support at large $x$.

We present a calculation of the pion quark momentum fraction, x , and its third Mellin moment x 2 . We also obtain, for the first time, x and x 2 for the kaon. We use an ensemble of two degenerate light, a strange and a charm quark (N f = 2 + 1 + 1) of maximally twisted mass fermions with clover improvement. The quark masses are chosen so that they reproduce a pion mass of about 260 MeV, and a kaon mass of 530 MeV. The lattice spacing of the ensemble is 0.093 fm and the lattice has a spatial extent of 3 fm. We analyze several values of the source-sink time separation within the range of 1.12 − 2.23 fm to study and eliminate excited-states contributions. The necessary renormalization functions are calculated non-perturbatively in the RI scheme, and are converted to the MS scheme at a scale of 2 GeV. The final values for the momentum fraction are x π u + = 0.261(3)stat (6)syst, x K u + = 0.246(2)stat(2)syst, and x K s + = 0.317(2)stat(1)syst. For the third Mellin moments we find x 2 π u + = 0.082 (21)stat (17)syst, x 2 K u + = 0.093(5)stat(3)syst, and x 2 K s + = 0.134(5)stat (2)syst. The reported systematic uncertainties are due to excited-state contamination. We also give the ratio x 2 / x which is an indication of how quickly the PDFs lose support at large x.

I. INTRODUCTION
The pion and kaon provide a good laboratory for studying QCD dynamics at hadronic scales. Moments of parton distribution functions (PDFs) are important quantities for the study of the internal structure of hadrons. They can directly be computed non-perturbatively within lattice QCD, using local operators, up to x 3 , and yield important insights that complement experimental programs that mostly measure PDFs for the nucleon. As the global initiatives to study PDFs in a variety of high-energy processes, such as deep-inelastic lepton scattering (DIS) and Drell-Yan in hadron-hadron collisions, at facilities such as Jefferson Lab, RHIC, Fermilab, and the LHC, intensify, providing theoretical insights on the moments has become very timely. In particular, studying the the pion and kaon moments will provide valuable information for the experimental program of the future Electron-Ion Collider [1,2].
PDFs provide a complementary picture of the structure of hadrons, as compared to electromagnetic form factors. However, unlike form factors PDFs are light-cone dominated quantities, and thus, cannot be calculated directly on a Euclidean lattice. A very promising approach to access the x-dependence of PDFs in lattice QCD has been proposed [3], and is developing in parallel with calculations of moments of PDFs. Significant progress needs to be made before reaching the precision and control of systematic uncertainties needed for a reliable comparison with experimental measurements. 1 The traditional approach to extract information on PDFs from lattice QCD is to evaluate Mellin moments of PDFs. Computing the Mellin moments of PDFs allows us to obtain results with more controlled systematics. They are interesting in their own right, as they are extracted from phenomenological analyses of experimental data, enabling direct comparison.
Pion and kaon structure is relevant to a number of important questions, such as, how hadron masses are generated, the dynamics of chiral symmetry breaking, and the role of pions in nucleon-nucleon interactions. The constrast between the nucleon, and the pion and kaon, is crucial to understand the Standard Model mechanisms that produce hadron masses. For example, in the chiral limit the masses of the pion and kaon vanish, whereas, the nucleon still has a mass on the order of 1 GeV. As such, the trace anomaly must vanish in the pion/kaon in the chiral limit but is non-vanishing in the nucleon. The pion, the lightest hadronic state in the QCD spectrum, is relevant for chiral symmetry breaking involved in nucleon-nucleon interactions. Pion cloud for instance, can explain why there are more d than u anti-quarks in the proton sea [27][28][29][30]. Comparisons between pion and kaon structure can reveal interesting aspects of QCD dynamics. For example, a model calculation [15] suggests that the strange contribution to the kaon form factors drops faster with increasing momentum transfer, compared to the up quark form factor, which has been interpreted as a consequence of confinement.
The rest of the paper is organized as follows: In Sec. II we present the theoretical setup of the calculation and the appropriate decomposition to obtain x and x 2 for mesons. In Sec. III A we provide the details of the lattice formulation, the parameters of the ensemble employed, as well as the calculation of the correlation functions needed in this work. The analysis for the extraction of reliable estimates for the non-perturbative renormalization functions is described in Sec. III B, as well as the details of the non-perturbative prescription. The various analyses on the two-point correlation functions for extracting the pion and kaon masses corresponding to the ensemble under study are presented in Sec. IV. In the same section we also include a thorough investigation of excited states for both x and x 2 , as well as an alternative kinematic setup for extracting x . The final results along with comparison with other studies and phenomenology are presented in Sec. V. In Sec. VI we summarize our work and conclude.

II. THEORETICAL SETUP
To calculate x , we compute the meson matrix elements M (p )|O {µν} |M (p) for p = p, where O {µν} is the one-derivative vector operator defined as where the notation {· · · } means symmetrization and traceless, with [31], and ψ corresponds to the up or down quark for the pion, and the up or strange quark for the kaon. {...} indicates symmetrization over the indices, in this case µ and ν, as well as subtraction of the trace, to avoid mixing with other operators [32][33][34].
The higher moment x 2 is accessed using a fermion operator with two covariant derivatives, O µνρ = ψγ µ D ν D ρ ψ. To avoid any mixing, we choose the indices µ, ν, ρ to be different [33,35,36]. Therefore, only a symmetrization over these indices is needed, that is, The meson matrix elements decompose into two generalized form factors, A 20 and B 20 for the one-derivative vector operator, and A 30 and B 30 for the two-derivative operator. The kinematic coefficients in Euclidean space are given by: In the above decompositions, P is the average of the initial and final momenta of the meson, P = (p+p )/2, and ∆ their difference, ∆ = p − p. Q 2 is the momentum transferred squared. C is a kinematic factor related to the normalization of the two-point functions. Therefore, C depends on the frame employed and the momentum transferred. Based on our conventions, we obtain C = for a general frame. m M is the mass of meson M and E(p)= m 2 M + p 2 is the energy at momentum p. Therefore, C does not depend on the spatial directions of the momentum, only on p 2 . Note also, that the generalized form factors A i0 and B i0 are functions of the momentum transferred squared, and are independent of the kinematic setup.
The quantities of interest are obtained from the forward limit of the matrix elements, leading to x ≡ A 20 (0), x 2 ≡ A 30 (0). The decomposition of Eq. (3) takes a simple form for mesons at rest, and in fact, at Q 2 = 0 there is only one matrix element contributing, which has µ = ν = 4 where we denote by index 4 the temporal direction. In such a simplified case, Eq. (3) becomes The index "M " indicates the meson of interest. The kinematic coefficient of A 30 becomes zero in the rest frame ( p = p = 0) and in the forward limit, unless all the indices of the operator are temporal. This is is not an optimal option, as O 444 suffers from mixing with lower-dimension operators [32][33][34], making the extraction of x 2 unreliable. In fact, to eliminate mixing even with operators of equal dimension, all indices must be different from each other, which is the choice we employ in this work. Given these constraints, the only way to obtain x 2 is to work in a frame in which the meson is moving with some momentum p = p = 0 (boosted frame). In the forward limit, the momentum is the same at the source and the sink, p = p = (iE, p x , p y , p z ). For µ = ν = ρ = µ, all spatial components of the momentum must be non-zero to extract x 2 directly from lattice data, without the need of applying fits on matrix elements with finite momentum transfer. In the boosted frame, the matrix elements are related to their corresponding Mellin moments via: which include all kinematic factors and normalizations. In Eq. (7) we take one of the indices to be temporal, which simplifies the kinematic factor of x 2 .

A. Lattice Action
In this work we employ one ensemble [37] of N f = 2 + 1 + 1 twisted mass fermions with a clover term and the Iwasaki improved gluon action. The ensemble is generated by the Extended Twisted Mass Collaboration (ETMC). The fermionic part of the action is written in the physical basis as D = γ µ (∇ * µ + ∇ µ )/2, ∇ µ and ∇ * µ are the forward and backward lattice covariant derivatives, and µ q is the twisted quark mass [26]. W cr = −(a/2)∇ * µ ∇ µ + m cr and m cr is the bare untwisted mass tuned to its critical value, which gives automatic O(a) improvement [38], requiring no further improvements on the operator level. The last term is the clover term multiplied by the Sheikoleslami-Wohlert improvement coefficient c sw . Since we achieve O(a) improvement from the critical mass, c sw is used to reduce isospin symmetry breaking effects [39]. We take c sw = 1.74 for this ensemble. The parameters of the simulation are given in Table I For the interpolating fields, J M , of the mesons under study we take A useful consequence of the pseudoscalar structure of the pion, as well as the γ 5 -hermiticity relation of the twisted mass quark propagators is that we only need to calculate the up quark contribution to the pion three-point functions. The pion and kaon interpolating fields are smeared using Gaussian smearing at both source and sink. The smearing parameters are tuned separately for the pion and kaon. We use the same value of α G but varying the number of smear iterations N G for the light and strange quarks. An optimal choice for N G is based on the criterion that the root mean squared radius of the smeared source reproduces the experimental radius of the pion [40] for the light quarks and the experimental radius of the kaon [41] for the strange quarks. In this work we obtain (α G , N G ) = (0.2, 50) for the light quarks and (α G , N G ) = (0.2, 40) for the strange quark. APE smearing is applied on the gauge links that enter the gaussian smearing with parameters (α AP E , N AP E ) = (0.5, 50). We study the connected contribution to the matrix elements of O 44 and O µν4 , which is shown in Fig. 1. For the calculation of the three-point functions we use the fixed sink sequential inversion approach. The three-point correlation functions are calculated at zero momentum transfer, where t i , t, t s are the source, insertion and sink Euclidean times, respectively. The corresponding spatial coordinates of the source, current insertion and sink are x i , x, x s . Without loss of generality we will take the source to be at t i = 0, so that the source-sink separation t s − t i = t s . For a general insertion current O Γ =ūΓu ±dΓd, the three-point functions can be written in terms of the up and down parts Performing the Wick contractions for the pion, and applying the γ 5 -hermiticity, it can be shown that, for π + , the up and down parts are related by The plus / minus sign comes from the fact that a general γ-structure is either γ 5 -hermitian or anti-γ 5 -hermitian, i.e. γ 5 Γ † γ 5 = ±Γ. Both the one-derivative vector and two-derivative vector operators are γ 5 -hermitian, and therefore, In the results presented here, we focus on the u + contribution to the pion, where the q + ≡ q +q notation has been adopted. Note that, based on Eqs. (15) and (16), x π u + = x π d + and x 2 π u + = x 2 π d + , that is, x π u + +d + = 2 x π u + and x 2 π u + +d + = 2 x 2 π u + . This discussion is relevant to the comparison with phenomenological results presented in Sec. VI. We note that Eqs. 14, 15, 16 are only applicable for the pion case, whereas for the kaon due to different mass of the light and strange quarks such relations do not hold.
We analyze 122 configurations, separated by 20 trajectories to reduce auto-correlation effects. In the rest frame, we use 16 randomly chosen source positions on each configuration, giving a total statistics of 1952. In the boosted frame, we use 32 source position for a total statistics of 3904. For the calculation in the rest frame, we use six source-sink time separations, namely t s /a = 12, 14, 16, 18, 20, 24, corresponding to t s = 1.12 − 2.23 fm. This allows for a thorough investigation and elimination of possible contributions from excited states. Based on the analysis of the results in the rest frame we concluded that a subset of t s /a = 14, 16, 18 is sufficient for extracting the ground state matrix elements. Thus, we only use these three time separations for the computation of x and x 2 , in the boosted frame.
According to the decomposition of Eq. (6), x 2 can be obtained using momentum boost with at least two nonzero spatial components, with the lowest momentum being p i = 2π L (±1, ±1, 0) (12 combinations). In this work, we employ, for x 2 , momenta of the class p i 2 = 12π 2 L 2 , which corresponds to 8 combinations for the spatial components, With the same setup we also obtain x , for a qualitative comparison with the rest frame, and the scaling of the statistical uncertainties. The choice p i 2 = 12π 2 L 2 is optimal for two reasons: While it increases the statistical uncertainties as compared to momenta p i 2 = 8π 2 L 2 , the computational cost for the same number of configurations is reduced by 33% due to the smaller number of permutations. Also, the class p i 2 = 12π 2 L 2 allows one to access, with the same setup, other quantities, such as x 3 , as well as form factors and generalized form factors. These quantities will be presented in a follow-up work.

B. Renormalization
The renormalization of the bare matrix elements is multiplicative, and the renormalization functions are calculated non-perturbatively using the Rome-Southampton method (RI scheme) [42]. The estimates are converted to the MSscheme and evolved at a renormalization scale of µ =2 GeV. We refer to the renormalization function of O µµ and O µνρ (µ = ν = ρ = µ) as Z vD and Z vDD , respectively. The renormalization function in the RI scheme are determined by the conditions where p is the momentum of the vertex function, set to the RI renormalization scale, is the treelevel value of the fermion propagator (operator), and the trace is taken over spin and color. We use the momentum source method [43], which is successfully employed for twisted mass fermions [44][45][46]. This method achieves per mil accuracy even with a small number of configurations. In the results presented here we use 10 configurations. In order to reduce discretization effects we use momenta that have the same spatial components, that is: where L t (L s ) is the temporal (spatial) extent of the lattice. These momenta are chosen to have suppressed non- , which is based on empirical arguments [47]. We improve the non-perturbative estimates for Z vD by subtracting finite lattice effects using the procedure outlined in Refs. [46,48]. The latter are computed to one loop in perturbation theory and to all orders in the lattice spacing, O(g 2 a ∞ ). Such a procedure is not yet available for two-derivative operators. However, we partly improve Z vDD , by subtracting the O(g 2 a ∞ ) artifacts from Z q which enters the renormalization condition for Z vDD in Eq. (17). For a proper chiral extrapolation, we calculate the renormalization functions on several ensembles with all masses of quark flavors degenerate (N f = 4). We use five ensembles at different values for the pion mass, which are produced with the same β value as the one of the cA211.32 ensemble analysed for the matrix elements. The parameters of the N f = 4 ensembles are given in Table II. The chiral limit is taken using a quadratic fit with respect to the pion mass. For both Z vD and Z vDD , we find a negligible dependence on the pion mass, as can be seen in Table III for two representative renormalization scales ((aµ 0 ) 2 = 2, 4). On each N f = 4 ensemble we use 23 values of the initial RI scale µ 0 ranging from ((aµ) 2 ∈ [0.9 − 6.7]). Each value is converted and evolved to MS(2 GeV) using an intermediate Renormalization Group Invariant scheme defined in continuum perturbation theory. A linear fit with respect to (aµ 0 ) 2 is applied on the MS values to eliminate residual dependence on the initial scale µ 0 . Such a dependence may be present due to finite-a effects and/or truncation of the conversion factor (to three loops in perturbation theory).
In Fig. 2 we show Z vD and Z vDD in the RI and MS schemes as a function of the initial RI renormalization scale, µ 0 . Z MS O are given at µ = 2 GeV, and the purely non-perturbative data exhibit a small residual dependence on the initial scale µ 0 they were evolved from. A procedure of subtracting the finite-a effects to O(g 2 a ∞ ) is also applied on Z vD . For Z vDD the improvement is only applied to Z q . We find that for both cases, subtracted results have a much smaller slope than the non-substracted ones, demonstrating the effectiveness of the artifact-subtraction procedure.
We eliminate any residual (aµ 0 ) 2 dependence in each renormalization function by using the Ansatz

A. Effective Mass
One of the important ingredients in the determination of the Mellin moments is the mass (energy) of the meson in the rest (boosted) frame. This is needed because the ground-state energy enters in the decomposition of Eqs. (3) - (6). We implement two fits for extracting the ground state energy from the two-point correlation functions as described below. We exploit the symmetry properties and we symmetrize the correlator corresponding to t and T − t, for t ∈ [0, T /2], i.e., the value at t have been averaged with their corresponding value at T − t.
Plateau method: The first method relies on a single-state fit where the effective mass (energy) is fitted to a constant with respect to t. The fit is taken over a range of t where the effective mass (energy) becomes time independent (plateau region). We calculate the effective mass from the symmetrized two-point function according to We test several values for the lower value of t entering the plateau fit, t low /a ∈ [11 − 19] for the rest frame and t low /a ∈ [8 − 14] for the boosted frame, while the maximum value is fixed to t/a = 31.
Two-state fit: The second method is a two-state fit in which we include the first excited state in the fit Ansatz given by The amplitudes c 0 and c 1 , as well as the ground and first excited state energies E 0 , E 1 are fit parameters. An alternative procedure is to apply the two-state fit on m eff directly, by substituting Eq. (22) into Eq. (21). This way, one of the amplitudes cancels out and the fit consists of three free parameters. Here we employ both procedures to cross-check the consistency of the results. We note that extracting the amplitude c 0 is needed in order to calculate the ratios between the three-point and two-point functions as described in the next section (see, e.g., Eq. (32)). The two-state fit is taken over the range t ∈ [t low − 31a], with t low /a ∈ [1 − 4]. In Fig. 3 we show the pion and kaon mass as a function of the lowest value of t/a entering the fit. We note that for the kaon we use the so-called Osterwalder Seiler (OS) fermions [49] which avoids the mixing effects between strange and charm quarks. The value of the µ s = 0.022 which enters the strange quark propagator is fixed by the physical ratio m Ds /f Ds = 7.9 [50], where m Ds is the mass of the D s meson and f Ds it's decay constant. We compare the results extracted from the plateau and two-state fits of Eq. (22). We find that there is a very good agreement between the two methods, when t low /a ≥ 11 in the plateau fit.
We repeat a similar process of extracting the energy and amplitude of the ground state, using the data in the boosted frame. As explained above, we focus on meson momentum boost p i 2 = 12π 2 L 2 . To increase the accuracy of the results and improve the stability of the fit, we perform the various fits on the averaged two-point functions over the eight values of the momentum boost leading to the same p i 2 ( p i = 2π L (±1, ±1, ±1)). The results of the fit are shown in Fig. 4 0  The notation is the same as in Fig. 3.
The final values shown with purple and blue in Figs. 3 and 4 are selected based on the following criterion: We accept a plateau fit with t low when the lowest state mass in the rest (boosted) frame (energy) obtained using the plateau method, m plat (E 0,plat ) and the the two-state fit m 2st (E 0,2st )satisfy the conditions where δm plat is the statistical error on the value extracted from the plateau method. An additional constraint for the accepted fit, is χ 2 plat /d.o.f < 1. Our final values for the pion mass in the rest frame based on the above criteria are: plateau : am π = 0.1250(2), t low /a = 11 , 2−state : am π = 0.1251(2), t low /a = 2 (25) while in the boosted frame we obtain plateau : aE 0π = 0.361(4), t low /a = 8 , A similar procedure for the kaon leads to plateau : am K = 0.2507(3), t low /a = 11 , 2−state : am K = 0.2508(2), t low /a = 2 (29) in the rest frame, and to plateau : aE 0K = 0.423(1), t low /a = 8 , (30) 2−state : aE 0K = 0.423(1), t low /a = 2 (31) in the boosted frame.
In Fig. 5 we plot the effective mass in the rest frame calculated from Eq. (21). We also show the plateau fit value and the two-state fit on the correlator as chosen based on Eq. (23). We find full agreement between the two fits, for both the pion and the kaon.

B. Excited-states contamination in x
To extract the ground state contributions to x , one has to ensure suppression of excited states in the three-point functions. This is achieved at sufficiently large insertion (t/a 1) and sink times ((t s − t)/a 1), where the ground state of the hadron gives the dominant contribution to the three-point correlation functions. It is in this region that we need to extract the matrix elements in order to control excited-states contamination. We employ six values of t s in the rest frame, t s /a ∈ [12,24], which for mesons can be achieved with a reasonable computational cost. For the pion this is due to the fact that the statistical error for meson matrix elements in the rest frame remains the same with increasing t s . Similarly to the analysis of the two-point functions we use two different analysis methods on the three-point functions, in order to study the convergence to the ground state and the significance of excited-states contributions. We used this study in the rest frame, as a guidance for the t s values to be employed in the boosted frame. Conclusions from such a study are also useful in a follow-up work for other studies of pion and kaon matrix elements.
Plateau method: The first method is based on a constant fit applied to an appropriate ratio of three-point and two-point functions. We choose a convenient ratio so that the denominator contains the ground state obtained from the fit of Eq. (22) (instead of the actual two-point functions). This removes the t s dependence in the ratio, allowing plateau convergence with increasing t s : We perform a constant fit as a function of the time t of the operator insertion for each t s separately and in a region where mild t-dependence is observed. One then seeks convergence of the extracted plateau values as t s increases. In the limit of large time separations the ratio becomes time-independent, that is Combining Eq. (5) with Eq. (33), we can obtain x via In the above expression we include the renormalization function for the one-derivative operator, Z vD , and all kinematic and normalization factors.
Two-state method: In the second method of extracting x we take into account the contribution from the firstexcited state in the three-point correlation functions. A two-state fit may be performed via where θ(t) = 1 for t ≥ 0 and θ(t) = 0 for t < 0. Given the large number of parameters, in Eq. (35) we use m for the rest frame and E 0 for the boosted frame and E 1 for the excited states extracted from the two-state fit of Eq. (22). Therefore, the actual free parameters are the amplitudes A 00 , A 01 , and A 11 (A 01 = A 10 for zero momentum transfer). The desirable matrix element of the ground state is extracted via where c 0 is the coefficient obtained from Eq. (22). Eq. (36) leads to the following expression for the renormalized x In Table IV we collect the values for x for the pion and kaon extracted from different source-sink time separations, and the two-state fit using t s /a = 12 − 24. The results for the two-state fit using t s /a = 14 − 18 is also included. The latter choice is based on investigating the dependence of x on the fit range. We find that the excited-state fit is compatible with the values obtained from t s /a 18 for both the pion and the kaon. As expected in the rest frame, the statistical uncertainties remain constant with increase of the source-sink separation. We find that the plateau values have statistical errors of 2% or less.
The ratios of three-to two-point functions for each value of t s are shown in Fig. 6, for the up contribution to the pion and the up and strange contribution of the kaon. R denotes the ratio R 44 M of Eq. 32 multiplied by all kinematic factors and the renormalization function. We observe that the excited-states contamination is similar for both the pion and kaon. We find convergence on x for t s 18a. A comparison of the two-state fit using t s ∈ [12a − 24a] and the plateau values is shown in the right panels of Fig. 6. Since the ground-state contribution is established at t s ≥ 18a and the two-state fits using t s ∈ [12a − 24a] and t s ∈ [14a − 18a] yield the same values, we choose to limit our calculations to t s /a = 14, 16, 18 for the boosted frame. A comparison for x between the two frames is discussed in Sec. IV C.  6. Results for the ratio leading to x . We show the up part of the pion, the up and strange parts of the kaon, in the top, center and bottom panels, respectively. In the left column, the points are the ratios in Eq. (32) multiplied by the renormalization and all kinematic factors. We plot values for ts/a = 12 up to ts/a = 24. The blue, red, green, magenta, cyan and orange constant bands are the results obtained from Eq. (34). The purple band is the two-state fit value obtained from Eq. (37). In the right column we plot the plateau values, together with two-state fit results. The gray band is the function obtained from a two-state fit using ts ∈ [12a − 24a] and taking t = ts/2.

C. Alternative setup for x in the boosted frame
In the above discussion we have used the rest frame for the extraction of x , p f = 0, and in the forward limit we also have p i = 0. In this paragraph we explore an alternative setup, a boosted frame with p f = p i = 0. Note that employing such a frame is not necessary for x , as A 20 (0) has a nonzero kinematic coefficient in the rest frame. However, the study of x within the boosted frame is interesting because one can understand how the statistical errors increase with t s . Based on the conclusions from Sec. IV B on the analysis of excited states, we focus on t s /a = 14, 16, 18, as the computational cost for the same number of configurations is increased by a factor of 8 as compared to the rest frame. Since this calculation is part of a wider set of operators, the optimal class of momenta is p 2 i = 12π 2 L 2 . This corresponds to eight combinations for the spatial components, that is, p = 2π L (±1, ±1, ±1). In such a frame, the appropriate decomposition is given in Eq. (6), instead of Eq. (5).
In Fig. 7 we compare the ratios leading to x for both the pion (top panels) and kaon (center and bottom panels). The left, center and right columns correspond to t s = 14a, 16a, 18a, respectively. The ratios include all the kinematic factors, and thus can be compared to each other. As can be seen, the statistical uncertainties increase with t s in the boosted frame, which is expected. For both the pion and kaon we find agreement between the plateau values obtained from the two frames within the uncertainties. The ratio for the x for the pion and the kaon is shown in Fig. 8 for t s /a = 14, 16, 18, and also compared to the two-state results. We find that all plateau values are compatible with the results of the two-state, indicating that excited-states contamination are within the reported uncertainties, which are larger compared to the ones in the rest frame. In Table V we collect all the results obtained in the boosted frame. We find that the statistical uncertainties in x π u grow from 5% to 10%, with the increase of t s from 14a to 18a. The corresponding increase in x k u ( x k s ) is 2% to 3% (2% to 4%). We remind that the error in the rest frame is less or equal to 2%, and it is constant for all source-sink separations. The extraction of x 2 is more challenging than x for several reasons. Firstly, x 2 cannot be extracted in the rest frame due to a vanishing kinematic factor in Eq. (7). The introduction of momentum in the meson states increases the statistical noise, and in our case, the use of a rather large momentum ( p 2 i = 12π 2 L 2 ) worsens the signal even more. Secondly, the presence of two covariant derivatives in the operator contribute to the increase of the gauge noise. Thirdly, having three Dirac indices leads to a more complicated renormalization pattern, and, to completely avoid the mixing with operators of equal or lower dimension, the indices of the operator must be selected different from each other. Here we employ the operator O µν4 .
In Fig. 9, we show the ratio leading to x 2 for the pion and the kaon. We plot the data for the three values of t s considered, that is t s /a = 14, 16, 18, and compare with the two-state results. We find that all plateau values are compatible with the results of the two-state, indicating that excited-states contamination are mild compared to the errors on this quantity. The values obtained from the plateau and two-state fits are given in Table VI

V. FINAL RESULTS AND COMPARISON WITH OTHER STUDIES
In this section we discuss our final values for the quantities studied in this work. For x we give the results in the rest frame and using the two-state fits. This is because the statistical uncertainties are the same for all values of t s . For x 2 we use the results extracted from t s /a = 18, as the two-state fit may be driven by the most accurate data. and These results are in the MS scheme at a scale of 2 GeV. We use the notation q + ≡ q +q for the sum from quark and antiquark contributions. As already mentioned x π u + +d + = 2 x π u + and x 2 π u + +d + = 2 x 2 π u + . The number given in the first parenthesis is the statistical error obtained from a jackknife analysis. We also report a systematic error, given in second parenthesis, which is due to excited states contamination. This is computed as the difference between the value extracted using the two-state fit and the plateau method for t s = 24 for x and t s = 18 for x 2 . We also extract the ratio x 2 / x , for which we find using the results at t s /a = 18 for x and x 2 . The number in the first parenthesis is statistical, while in the second parenthesis is systematic due to excited states.
There is very limited experimental data on the kaon PDF, so it is interesting to contrast these moment results with expectations from model calculations. For our lattice results we find x K u + < x π u + < x K s + which is consistent with many phenomenological calculations, including the DSE results of Ref. [14]. This ordering in the momentum fractions is understood because the heavier s quark skews s K (x) to larger x which is compensated by a shift in u K (x) to smaller x. In the limit of equal quark masses these moments would be degenerate, therefore, we find flavor breaking effects of up to 20% in these moments. For the third Mellin moment we find x 2 K u + , x 2 π u + < x K s + , which is again consistent with expectations. Uncertainties on the third Mellin for the u quark in the pion and kaon do not allow an ordering of these moments, however, any deviation from the order found for the momentum fractions would be very interesting.
There are a number of calculations on the pion x [22,24,26,51,52], including results obtained directly at the physical point [25]. The pion third Mellin moment x 2 , on the other hand, is lesser known, and has been studied in Refs. [22,26,51] using different lattice formulations. Here we compare with lattice results on x π u + extracted at the same or similar value of the pion mass, that is 240 -270 MeV.
It is interesting to compare with phenomenological estimates, which can be found in Refs. [53] and [54]. The older analysis of Ref. [54] gives a value of x π u = 0.217 (10), in the MS at a scale at (5.2 GeV) 2 , while ours is (2 GeV) 2 . Converting to 2 GeV, their value becomes x π u = 0.361 (17). A more recent analysis is presented by the JAM Collaboration [53], on a large set of experimental data including Drell-Yan data and leading neutron electroproduction from HERA. They find x π valence = 0.480 (10), which is reasonably close to our value of 2 x π u + = 0.522 (13). The error in the parenthesis is the combined statistical and systematic uncertainties added in quadrature. The fact that our value is higher, by ∼ 4%, maybe attributed to the fact that our calculation is not at the physical point and the continuum limit is not taken. Both the chiral extrapolation and taking a → 0 will decrease this value as demonstrated in Ref. [26]. We also note that all lattice calculations to date consider only the connected contributions as done in this work. The disconnected contributions should be included for a final comparison with phenomenology. We summarize the results for x π u in Table VII. There are very limited calculations for x 2 within lattice QCD, and the one which we can directly compare with our results is Ref. [26]. They find x 2 π u + = 0.131(18)(24) and x 2 π u + = 0.132(40)(53) for ensembles A30.32 and B25.32, respectively. These estimates are compatible with our final value, within the large uncertainties of the aforementioned values. x 2 π was also calculated in Refs. [22,51] using a different operator, which has two Dirac indices the same. Such a choice is expected to lead to more complicated renormalization pattern due to mixing, which is not addressed in Refs. [22,51]. They obtain x 2 π u + = 0.128(9)(4) which is, however consistent with the value obtained in this work. Phenomenological estimates for x 2 π can be found in Ref. [53] where a value of x 2 π valence = 0.210(5) is reported, which is compatible with our value of 2 x 2 π u + = 0.164(54) within uncertainties. However, one needs to bear in mind that the phenomenological value does not include sea quark contributions unlike the lattice QCD calculation, where such sea quark effects are automatically included. For completeness, we also provide the results from Ref. [54], which correspond to a scale of (5.2 GeV) 2 . Their finding is x 2 π u = 0.087 (5). We convert this result to 2 GeV resulting to x 2 π u = 0.169 (10). This is compatible with our value.

VI. SUMMARY
We present a calculation of the second and third Mellin moments, x and x 2 for the pion and kaon. We use one N f = 2 + 1 + 1 ensemble reproducing a pion mass of 260 MeV and a kaon mass of 530 MeV. For x we employ both the rest and boosted frames, and we find full agreement between the two. However, the statistical uncertainties for the boosted frame are larger, as can be seen in Fig. 7. To extract x 2 one requires a boosted frame due to kinematical factors. The selected meson momentum boost has all spatial components nonzero, and gives p i 2 = 12π 2 L 2 (0.52 GeV 2 ). We renormalize all matrix elements with non-perturbative renormalization with cut-off subtraction that utilizes lattice QCD perturbation theory.
We perform a thorough investigation of excited states using the three-point function that determines x . For this investigation we use the rest frame and calculated the matrix elements for six values of the source-sink time separation ranging from 1.12 fm to 2.23 fm. The computational cost for this study is within reach, as the statistical error does not increase with t s in the rest frame for the pion and only increases mildly for the kaon. We analyze the data using one-state and two-state fits. We find that excited states are suppressed for t s > 1.6 fm, that is, t s = 18a or higher. Another important conclusion from the analysis is that the two-state fits using only t s /a = 14, 16, 18 are compatible with the two-state fits obtained including the larger t s values. This is crucial, as in the case of the boosted frame, the statistical errors increase with t s , as illustrated in e.g, Fig. 9, limiting how large t s can be. Thus, for the boosted frame we perform the computation for t s /a = 14, 16, 18, where the consistency of the results extracted between oneand two-state fits demonstrates that excited states are correctly accounted for.
The results for the pion are given in Eq. (38) and Eq. (41), in the MS scheme at a renormalization scale of 2 GeV. Our results agree very well with the lattice QCD analysis of Ref. [26]. It is important to emphasize that the in-depth study and elimination of excited-states in our analysis has reduced the extracted values bringing them closer to those determined from phenomenology. For example, our lattice data for source-sink time separations below 1.6 fm give a value that is 10% − 20% higher than the final value extracted when the larger separations are used (see Table IV).
Our final results for x K u,s and x 2 K u,s are given in Eqs. (39) - (40) and Eqs. (42) - (43), respectively in the MS scheme at a scale of 2 GeV. Currently, there are no other lattice data for these quantities, nor global fits on experimental data. Therefore, the results on the kaon presented in this work provide a first prediction. Taking in to account that for the ensemble employed in this work the kaon mass is about 530 MeV, i.e. only ∼ 7% heavier than its physical value, means that the values for x K u,s and x 2 K u,s can be considered as a good approximation of their physical counterparts.
In the near future, we will consider calculation of x 3 for both the pion and kaon. Another direction is the form factors and generalized form factors, which require off-forward matrix elements. We intend to test the momentum smearing method [55] for the boosted frame, which has been proven to increase the overlap with the ground state, decreasing significantly the statistical noise. This has already been implemented by members of our group in nucleon matrix elements of non-local operators [56][57][58][59][60][61], and can easily be applied for meson matrix elements.