$\langle x\rangle$ and $\langle x^2\rangle$ of the pion PDF from Lattice QCD with $N_f=2+1+1$ dynamical quark flavours

Using Nf=2+1+1 lattice QCD, we determine the fermionic connected contributions to the second and third moment of the pion PDF. Based on gauge configurations from the European Twisted Mass Collaboration, chiral and continuum extrapolations are performed using pion masses in the range of 230 to 500 MeV and three values of the lattice spacing. Finite volume effects are investigated using different volumes. In order to avoid mixing under renormalisation for the third moment, we use an operator with two non-zero spatial components of momentum. Momenta are injected using twisted boundary conditions. Our final values read $\langle x\rangle=0.2075(106)$ and $\langle x^2\rangle=0.163(33)$, determined at 2 GeV in the $\overline{MS}$-scheme and with systematic and statistical uncertainties summend in quadrature.

Finite volume effects are investigated using different volumes. In order to avoid mixing under renormalisation for the third moment, we use an operator with two non-zero spatial components of momentum. Momenta are injected using twisted boundary conditions. Our final values read x phys R = 0.2075(106) and x 2 phys R = 0.163 (33), determined at 2 GeV in the MS-scheme and with systematic and statistical uncertainties summend in quadrature.

I. INTRODUCTION
In Quantum Chromodynamics (QCD), the pion represents the Goldstone boson of spontaneously broken chiral symmetry and is the lightest hadronic state in the spectrum. As such it is of deep importance both for the long range part of the nucleon-nucleon interaction and for the inner structure of the nucleon. In the latter case, it is now widely recognized that the pion is responsible for most, arXiv:1810.09743v1 [hep-lat] 23 Oct 2018 to extract moments beyond the leading one. This progress has led to investigations of interesting physics questions, such as the momentum and spin decomposition of the nucleon in terms of their quark and gluon contents. Within our ETM collaboration, this has been accomplished by lattice QCD simulations directly at the physical point [27], where both momentum and spin sum rules have been verified without imposing any constraints. For the pion, however, the situation is much less satisfactory. Earlier studies have computed the first three moments either in the quenched approximation [28][29][30] or for connected insertions only [31][32][33][34], all of them using simulations with quark masses away from their physical value. Only a few results for realistic QCD simulations appear in the literature, that is, Ref. [33] at about 150 MeV for the pion mass, and a determination directly at the physical point with N f = 2 in Ref. [34] for the lowest moment. Given the importance of the pion for ongoing and planned experiments, further study of the pion structure is imperative.
For the extraction of reliable estimates, systematic uncertainties such as discretisation and volume effects must be properly addressed and quantified.
In this work we present a multi-component effort in the aforementioned direction, with a variety of improvements compared to the studies available in the literature, in terms of the ensembles employed and level of control over systematic uncertainties in the computed moments. We calculate the light quark connected contributions of the second and third moment -x and x 2 -of the pion using lattice QCD simulations that include degenerate light as well as strange and charm quarks in the sea (N f =2+1+1). We use several ensembles produced by the European Twisted Mass Collaboration (ETMC) corresponding to three values of the lattice spacing, which allow us to study discretisation effects. These ensembles have pion mass values that range between 230 MeV and 500 MeV, which are combined in a chiral extrapolation to obtain the value at the physical pion mass. Different volumes are used to investigate finite size effects and excited state contaminations. In addition, a way around possible mixing for x 2 is the choice of an operator that is free from mixing under renormalisation. A first preliminary account of this work can be found in Ref. [35].
The remainder of the paper is organized as follows: In Section II we discuss the technical aspects of the lattice calculation, while in Section III we discuss the method used for the determination of the required renormalisation functions in the RI'-MOM scheme and the conversion to the MSscheme. The main results of this work are presented in Section IV, followed by a discussion and a summary in Section V. Technical details related to renormalisation can be found in Appendix A, while correlation coefficients of fit parameters are collected in Appendix B. Tabelle I: The N f = 2 + 1 + 1 ensembles used in this investigation, where the notation of Ref. [38] is used for labeling. In addition to the bare parameters and lattice time and spatial extents T /a and L/a, respectively, we give the number of configurations N conf on which the moments have been estimated.

II. LATTICE DETAILS
The calculation presented in this paper is based on gauge configurations generated by ETMC with N f = 2 + 1 + 1 dynamical quark flavors at three values of the lattice spacing. Details for the configuration generation and analyses for basic quantities can be found in Refs. [36][37][38]. The ensembles were generated using the Iwasaki gauge action [39] and the Wilson twisted mass fermion action at maximal twist [40][41][42]. Working at maximal twist guarantees O(a) improvement for most physical quantities [40], and in particular for the quantities considered here.
The bare parameters of the ensembles used here are summarized in Table I. µ is the bare light quark mass directly proportional to the renormalised light quark mass. µ σ and µ δ parametrize the strange and charm quark masses [36,41]. For the subset of configurations we used from each ensemble we have computed the autocorrelation times for the relevant quantities to verify their statistical independence. The error analysis is performed using the stationary blocked bootstrap procedure [43] with 1500 bootstrap samples.
In general, the computation of the moments requires the computation of three-point functions of the form with operator O inserted at Euclidean time t. We fix here the time difference between the two pions to T /2, which is not necessary, but convenient. The operators for the two moments will be detailed below. The particular choice of operators is motivated by their transformation properties under the symmetries of the lattice as well as the requirement of minimal mixing with lower-dimensional operators under renormalisation, see Refs. [25,44,45] for details. The interpolating operators for the pions read with the momentum p ∝ θ − θ realized via twisted boundary conditions, see below. τ i , i = 1, 2, 3 are the Pauli matrices acting in flavor space and ψ = (u, d) t is the light quark field.
A. The second moment x A convenient operator in Euclidean space-time for the calculation of the second moment x is Here, ↔ Dµ= 1 2 ( µ + * µ ) is the symmetric, gauge covariant lattice derivative with µ ( * µ ) being the usual gauge covariant forward (backward) derivative on the lattice. The above operator has the advantage that x is extracted without the need for an external momentum, because external momentum in general increases the noise. For the first use of this operator with Wilson twisted mass fermions we refer to Refs. [30,46].
The bare moment x bare is related to the matrix element of the operator O 44 as follows where p = (p 0 , p) is the four momentum of the pions. With pions at rest one obtains with M π the mass of the pion. The matrix element π(0)|O 44 |π(0) between two pions at rest is calculated from the ratio π(0) π(T /2) of the three point function over the two point function In Eq. 6 a factor of 2M π relates the lattice and continuum matrix elements of O 44 between pion states and a further factor of 2 relates the ratio of correlation functions to the value of the matrix element. This leads to There are two contributions in the Wick contractions of C 44 : the first is extracted when the current couples to the quarks of the pion directly (connected diagram), while the second is obtained from the so-called quark loop (disconnected diagram) in which the current interacts with the pion via gluon exchange. Both are visualized in Fig. 1. The disconnected contribution is ignored in our calculation, assuming that it is small, which is indeed the case for the nucleon [47]. A computation of the disconnected contributions is, however, planned for the near future.
B. The third moment x 2 In order to avoid mixing under renormalisation with lower dimensional operators [44], we use for the third moment the operator (in Euclidean space-time) which is related to x 2 via π(p)|O 012 |π(p) = −2(p 0 p 1 p 2 ) x 2 bare .
In contrast to x , non-zero momentum is needed to extract x 2 , due to the presence of the kinematic factor (p 0 p 1 p 2 ) multiplying the quantity of interest. We use twisted boundary conditions here to inject momentum, see below. As in the case of the second moment, the matrix element of the third moment is related to a ratio of three point to two point function which leads to For details on the implementation of O 012 we refer to Ref. [44], noting that while we employ the convention given therein for the discretisation of terms involving to the data for C π (t) for sufficiently large Euclidean times. In twisted mass lattice QCD at maximal twist the pion decay constant is directly related to the amplitude A via without the need for renormalisation [49].
We note that M π and f π are affected significantly by finite size effects [48]. Therefore, we use the corrections computed in Ref. [48], which are summarized in table II for all the ensembles used in this work.
In Eqs. 9 and 13 one needs to divide by the two-point function at T /2. We explore two possibilities to perform this division: the first one is to use the data of C π for t = T /2. The second one is to first fit Eq. 14 to the data for C π in the region, where the ground state dominates and then use the best fit parameters to reconstruct C π (T /2). The latter procedure can help to average out fluctuations.
However, the differences between the two procedures are always well below the statistical uncertainty both for x and x 2 . As our method of choice, we henceforth use the second of the two methods. Tabelle II: Per ensemble we give the θ-values for spatial twisted boundary conditions, pion masses and decay constants determined from C π as well as finite volume correction factors from Ref. [48].

D. Stochastic Evaluation
The above two and three point correlators are evaluated by using a stochastic time slice source (Z(2)-noise in both real and imaginary part) [50][51][52] for all color, spin and spatial indices. i.e., the for X, where the Z(2) random source ξ( z) a α satisfies the random average condition The generalized propagator [53] for Σ. This approach was first applied for x of the pion in Ref. [32] and we used it recently in a computation of the pion vector form factor [35,54]. To further improve the signal, we use N src = 5 sources per gauge configuration and average. The source time slices are chosen uniformly random

E. Twisted Boundary Conditions
In order to realize non-zero momentum of arbitrary values for the pions as needed for x 2 , we make use of so-called twisted boundary conditions [55][56][57]. Enforcing the spatial boundary conditions ψ(x+ e i L) = e 2πiθ i ψ(x) on the quark fields changes the momentum quantization condition in finite volume to p i = 2πθ i L + 2πn i L . In time direction we chose θ 0 = 1/2 to obtain anti periodic boundary conditions in time. We chose the θ in the spatial directions non-zero to obtain non-zero momentum for the pions.
For the two quarks in the pion, we always chose zero twist angle for one of the quarks and non-zero θ for the other one. The pion three-momentum p is then given by (n i = 0) We recall that for the computation of x 2 two non-zero spatial components of the pion momentum are needed when O 012 is used, see Eq. 13. We chose the two non-zero elements of θ equal, i.e. for instance θ = (θ, θ, 0). The corresponding values for θ for each ensemble are compiled in table II.
We always perform the computation for x 2 for both ± p and average. The such obtained result is The main reason for using twisted boundary conditions is the fact that noise in the three point and two point functions increases significantly with increasing modulus of the injected momentum.
With twisted boundary conditions we are able to chose the momentum as small as possible. However, we remark that twisted boundary conditions induce additional finite volume effects, which might influence our results [58]. As will be discussed later, we do not see such effects in x 2 within statistical uncertainties.

F. Chiral Extrapolations
In Ref. [59] the pion mass dependence of pion moments has been computed in leading order (LO) chiral perturbation theory (ChPT). The functional form for x reads with low energy constants (LECs) c 0 and c 1 . For the third moment it reads where we denote the corresponding LECs with b 0 and b 1 . We chose the renormalisation scale conventionally µ R = f π . In contrast to Ref. [59], we have expressed the two moments as a function of M π /f π , which has the big advantage of fully dimensionless expressions. In principle one should then use f phys π , i.e. the physical value of the decay constant. However, we use here f π as estimated for each ensemble, because scale setting is required only to estimate the moments at the physical point. Since f π is a constant in leading order ChPT, this procedure is consistent to the order of ChPT we are working here. Unfortunately, the next-to-leading-order expressions for the moments are not known. Still, in contrast to the case of nucleons, in the pion sector ChPT works well such that we expect already the lowest order to provide a reliable tool for our set of pion masses. In order to account also for lattice spacing artifacts we add terms c a a 2 /r 2 0 and b a a 2 /r 2 0 to the expressions for the second and the third moment, respectively.

III. RENORMALISATION FUNCTIONS
A renormalisation factor (Z-factor) must be applied to the bare matrix elements of the operators defined in Eq. (3) and Eq. (10), in order to obtain physical quantities. More precisely, the bare and the renormalised moments are related as follows In particular, the renormalisation procedure eliminates divergences with respect to the lattice regulator, and allows the continuum limit to be taken. In this Section we present the methodology and results for the renormalisation functions, which are finally converted to the MS-scheme at a scale µ=2 GeV. We employ the Rome-Southampton method (RI scheme) [60] to compute the Z-factors non-perturbatively determined by the conditions The momentum p is set to the RI renormalisation scale, µ 0 , S Born (Γ Born ) is the tree-level value of the fermion propagator (operator), and the trace is taken over spin and color indices.
We obtain the Z-factors using several ensembles at different values of the pion mass, so that the chiral limit can be safely taken. We employ the momentum source method introduced in Ref. [63] and used in Ref. [64], which leads to a high statistical accuracy with a small number of configurations. For the Z-factors presented in this work we use between 10 to 50 configurations depending on the ensemble under study. To reduce discretisation 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, and we restrict the momenta up to (a p) 2 ∼7. A useful constraint for the chosen spatial momenta is i p 4 i /( i p 2 i ) 2 <0.3 which ensures reduced discretisation effects. This is based on empirical arguments [65], as this ratio appears to suppress O(a 2 ) terms in the perturbative expressions for the Greens functions. The procedure we follow in this work is the same as our previous work in non-perturbative renormalisation, and thus, we refer the interested reader to Refs. [64,66,67] for technical details. It is worth mentioning that in the renormalisation of the one-derivative operator we also employ improvements by subtracting lattice artifacts [64]. The latter are computed to one-loop in perturbation theory and to all orders in the lattice spacing, O(g 2 a ∞ ). These artifacts are present in the non-perturbative vertex functions of the fermion propagator and fermion operators under study. Such an improvement is not yet available for the two-derivative operator, but finite a effects are partly removed upon the (ap) 2 →0 extrapolation. In this Section we focus on the results for Z vDD , which are presented for the first time, while results on Z vD have been extracted within the work of Ref. [64]. To extract the renormalisation functions in the chiral limit we perform an extrapolation using a quadratic fit with respect to the pion mass of the ensemble, that is, a RI (µ 0 ) + b RI (µ 0 ) · M 2 π where a and b depend on the scheme and scale. In addition, these parameters depend on the coupling constant, and separate fits are performed at each value of β. We find that the renormalisation functions under study have a very mild dependence on the pion mass, which becomes slightly larger for (a p) 2 <1. However, these points do not participate in the fit (ap) 2 →0 for the final estimates.
Allowing a slope, b =0, and performing a linear extrapolation with respect to M 2 π the data yield a slope that is compatible with zero within the small uncertainties. This is demonstrated in Fig. 2 for Z RI vDD for β=1.95, plotted as a function of the initial scale (ap) 2 . For clarity we only show two values of the twisted mass a µ sea , while the statistical errors are too small to be visible. The corresponding plot for Z vD is shown in Ref. [64]. In order to compare lattice values to experimental results one must convert to a universal renormalisation scheme and use a reference scale µ. Typically one chooses the MS-scheme at µ=2 GeV.

The conversion from RI to MS scheme uses the intermediate Renormalization Group Invariant
(RGI) scheme, which is scale independent and thus, The conversion factor can be extracted from the above relation .
The quantity ∆Z S O (µ 0 ) is expressed in terms of the β-function and the anomalous dimension γ S O ≡ γ S of the operator with all necessary ingredients defined in Appendix A. We employ a 3-loop approximation, for which ∆Z S O (µ 0 ) takes a simpler form [64]. In Fig. 3 we present representative results for Z vDD (at β=2.10) in the RI (Z RI vDD (µ 0 )) and MS (Z MS vDD (2 GeV)) schemes as a function of the initial RI renormalisation scale, µ 0 =p. Note that Z MS vDD has been evolved to 2 GeV, but there is residual dependence on the initial scale. This dependence is removed by extrapolating to zero, using the Ansatz where Z O corresponds to our final value on the renormalisation functions for the operator O. For each value of β we consider momenta 6 ≥ (a p) 2 ≥ 2 for which perturbation theory is trustworthy and lattice artifacts are still under control.
In Table IV we report our chirally extrapolated values for Z  In Ref. [68] Z vD has been determined on the same gauge configurations using a different method.
The authors find generally lower values for Z vD , which are within the quoted systematic uncertainties compatible with what we quote in Table IV. In order to be consistent in our treatment of x and x 2 , we stick here to the values compiled in Table IV.

IV. RESULTS
In the left panel of Figure 4 we show the bare three-point function C 44 (t) for ensemble D30.48.
The plot demonstrates the quality of the data we are able to obtain for C 44 (t). In Ref. [69] it was found that, in the quenched approximation, with Schrödinger functional boundary conditions and clover improved Wilson fermions, finite size effects to x are quite sizable. This persists even at values of M π · L where finite size effects in M π are no longer visible. The authors measured these effects to be of about 5% at values of M π · L ≈ 4. In this work we use two ensembles, A40.24 with M π · L = 3.5 and A40.32 with M π · L = 4.5, which differ only in the volume, and can be used for investigation of finite size effects. In the right panel of Figure 4 we show a comparison of x (t) between A40.24 and A40.32. We find that the values of the bare x (t) in the plateau regions for A40.24 and A40.32 agree within error bars. This indicates that in our lattice discretisation for the given values of M π · L, finite size effects play a minor role, if any. This is in agreement with the finding in the quenched approximation [46].
In Figure 5  For determining an estimate of x and x 2 we perform plateau fits to the (anti-) symmetrized data for x (t) and x 2 (t). Following the ideas put forward in Ref. [70], we perform such fits for many different fit ranges. The estimates for M π and f π are obtained by fitting all possible fit ranges with at least 6 consecutive time slices. Each of these fits obtains a weight according to Here, p i is the p-value of the corresponding fit and ∆ i is the statistical error on M π or f π determined from the bootstrap procedure for this fit range. This procedure is repeated for each bootstrap replica. In addition, a systematic uncertainty from the fit range choice can be specified from the 68% confidence interval of the weighted distribution.
The estimates for the moments are then obtained in a very similar manner, just that also M π is needed. Thus, we combine all possible fit ranges with at least six consecutive time slices for C π with all possible fit ranges with at least six consecutive time slices to the corresponding three point function. The weight for a moment with a specific fit range combination is obtained by multiplying the corresponding weights of the fit to C π and the fit to the three point function.  (14)(12) Tabelle V: The results for the renormalised x R and x 2 R for the N f = 2 + 1 + 1 ensembles used in this investigation. x R and x 2 R are given at 2 GeV in the MS-scheme. In addition we give the values of M π · L.  These results for the renormalised second and third moments of the pion are shown in Figure 6 in the left and right panel, respectively. They are plotted as a function of (M π /f π ) 2 with statistical errors only. Abbildung 6: x R of the pion for N f = 2 + 1 + 1 (a) and x 2 R (b) both renormalised atμ = 2 GeV in the MS scheme as functions of (M π /f π ) 2 . Dashed colored lines represent the best fit functions Eq. 30 at the three lattice spacing values, respectively. The solid black line represents the continuum curve. The black triangles represent the estimate at the physical point in the continuum limit. The error bars represent the statistical uncertainty.

A. Chiral and Continuum Extrapolations
The ChPT expressions Eq. 19 and 20 plus terms proportional to (a/r 0 ) 2 for x and x 2 read We perform fits of these functional forms to all the data of the second and third moment separately.
For these fits we have the data for the bare x ( x 2 ) and the estimates for Z vD (Z vDD ) and r 0 /a. To properly account for the uncertainties in the renormalisation functions and the Sommer parameter r 0 /a, we use the augmented χ 2 function as follows Here, Z and ∆Z denote the relevant renormalisation factor and its statistical uncertainty either for the second or the third moment. P Z and P r 0 are additional fit parameters per β-value. The usual χ 2 function entering χ 2 aug reads Here i(β) index the data points for the corresponding β-value, y i is the bare data for x ( x 2 ) and x i the data for (M π /f π ) 2 . With {P } we label the set of fit parameters {c 0 , c 1 , c a , P r 0 , P Z vD } ({b 0 , b 1 , b a , P r 0 , P Z vDD }) and with g the corresponding ChPT expression. The equation for the χ 2 function above is written for the uncorrelated case, because all data points stem from independent ensembles and r 0 /a and the renormalisation constants from independent analyses. Errors of fit parameters are again computed using the bootstrap procedure by performing a fit on every bootstrap replica.
In principle one could also include the error on (M π /f π ) 2 in the fit. However, these errors are so small compared to the ones for the moments that they do not alter the fit results. We also do not include systematic uncertainties in the fit, because they lack a statistical interpretation and would increase all error bars more or less uniformly.
The p-value of the fit equals 0.61 with χ 2 aug /dof = 8.2/10. Thus, the fit is acceptable and the continuum value of x at the physical point -defined via M π /f π = 1.0337 -reads x phys R = 0.2075 (53) .
The best fit curves for the three lattice spacings are included as dashed lines in the left panel of Figure 6. The continuum curve is the black solid line with the estimate of the second moment at the physical point indicated with the black triangle.
If we had used the values for Z vD from Ref. [68] instead of the ones we quote in Table IV, we had obtained an equally good fit with result x phys R = 0.189 (16). This value is compatible with the value above, but with larger uncertainties. We include half of the difference as systematic error in our final value.
For the third moment, the best fit parameters read b 0 = 0.16(2) , b 1 = 0.005(2) , b a = −1.6(7) With a p-value of 0.89 (χ 2 aug /dof = 5/10) the continuum estimate at the physical point reads Like for x R , the corresponding curves are shown in the right panel of Figure 6 in addition to the data. Again, the black triangle represents the estimate of the third moment at the physical point in the continuum limit.
For both the second and the third moment the fit parameters for the renormalisation factors and for r 0 /a agree very well with the input data. All best fit parameters, their uncertainties and correlations are compiled in Appendix B.

V. DISCUSSION
In this work we demonstrate the feasibility of the lattice calculation for the second and third moments of the pion PDF. Despite the challenges present in calculations of higher moments, we find sufficiently long plateau regions for the bare matrix elements for all ensembles used here, with the dependence on the fit-range of the order of the statistical uncertainties. For the second moment, where our bare values are precise to the few percent level, we observe a sizable dependence on M 2 π and significant lattice artifacts, cf. the left panel of Figure 6. From the value of x R of ensemble D15.48, which is the smallest pion mass value closest to the continuum limit, there is still a 10% difference to the continuum value at the physical point.
The statistical errors for x 2 are significantly larger than for x , since two derivatives and two non-zero spatial components of momentum are required. Therefore, pion mass and lattice spacing dependence are both not significant: all the data could be fitted to a constant in M π /f π with a result similar to the one we quote above. For both moments, finite size effects turn out to be not relevant, which is in agreement with the finding of Ref. [34], where the twisted mass formulation was used as well.
Our results for the second moment can be compared to other lattice computations, including our recent work using N f = 2 simulations directly at the physical point, however, without extrapolation to the continuum limit [34]. The value found in Ref. [34] also neglecting disconnected contributions at the physical point reads x R = 0.214(15)( +12 −9 ). It is fully compatible with the result we find here. In Refs. [31,71] a value of x R = 0.271(2)(10) atμ = 2 GeV in the MS scheme is quoted for N f = 2 flavor QCD also neglecting disconnected diagrams, which is significantly larger than our value. In these two references almost no lattice artifacts appear to be visible, in contrast to our findings. In the work of Bali et al. [33] a significantly lower value is reported, using a single ensemble at near physical pion mass value.
It is not so easy to identify a reason for the differences we observe. It seems the number of flavors is not so important, because our result with N f = 2 + 1 + 1 quark flavors is fully compatible with the N f = 2 result at the physical point. Even though the latter computation is at a single lattice spacing only, lattice spacing effects seem to be small with this action [72]. Thus, differences are likely to come from the chosen lattice discretisation leading to different lattice artifacts and finite size effects. This clearly demands furher careful investigations of systematic uncertainties in the future.
Refs. [31,71]  On the other hand, in practice the fermionic connected and eventually also the disconnected contributions can be determined. It is then very appealing to identify the part coming from the disconnected contributions as purely sea moments, see fig. 1. This allows one to make contact to the phenomenological point of view, where typically the following sum rule is used, here for x where the v, s, G denote the valence quark, sea quark, and gluon contributions, respectively. On the other hand, for a lattice calculation one would write: x conn where conn (disc) stands for a lattice computation performed with only fermionic connected (disconnected) contributions to the corresponding three-point function taken into account. The sum in q is over all active quark flavors. As defined in Eq. 3, the quantity calculated in this work is the total connected only contribution: Still, from Eq. 39 it is clear that x phys R (μ 2 ) cannot be the valence contribution of Eq. 37, because the connected contributions also receive contributions from so-called Z diagrams, which are counted as sea quark distributions in Eq. 37. Nevertheless, since the following equality must hold: we may, keeping the caveat discussed above in mind, compare x phys R (μ 2 ) with phenomenology if we understand the quantity computed here as an upper limit for x v (μ 2 ).
Phenomenological results for average x and x 2 are provided in Refs. [73] and [74]. Below we compare to the more recent results from Ref. [73], which are based on a larger set of experimental data, where they find  [77] to measure the pion structure functions will be instrumental to settle this matter, having also an impact in the decomposition of the pion momentum sum rule. Our value for x 2 is larger than 2 x 2 v , but its also has a large error bar.
Finally, we note that the relative share of connected to disconnected contribution to the total x may well depend on the pion mass.

VI. SUMMARY
In this paper we have presented results for the second and third moment of the pion PDF computed in N f = 2 + 1 + 1 lattice QCD. While we still neglect fermionic disconnected diagrams for both moments, we have thoroughly investigated the extrapolations to the physical point and to the continuum. This was possible due to ETMC ensembles spanning three values of the lattice spacing and pion masses ranging from 270 to 500 MeV. For x and x 2 we use operators which avoid any mixing under renormalisation. By studying two ensembles with all identical parameters but the lattice size, we can exclude finite volume effects significantly larger than our statistical uncertainties.
For the computation of x 2 non-zero spatial momenta are required which we inject using twisted boundary conditions. These allow us to chose the momenta optimally for the signal to noise ratio in the corresponding three-point function. Still, our results for x 2 have significantly larger statistical uncertainties than the ones for x , which is of course also due to the second derivative needed for It turns out that the choice of fit-ranges represents a major systematic uncertainty in the calculation of the moments. We approach this uncertainty by performing many fits and include them all weighted appropriately in the final estimates. From the weighted distribution a systematic error can be estimated which is typically of the order of the statistical error. The only exception is our ensemble at the smallest lattice spacing and pion mass value, where the systematic errors prevent us from obtaining a significant result. In summary we obtain determined at 2 GeV in the MS-scheme. In the bare matrix elements we find on average a 1% systematic error on x and a 15% systematic error on x 2 , which we have added to the final results in order to reflect the systematic uncertainty coming from the fit range choice. In x we add the systematic uncertainty from using the Z-factors determined in Ref. [68] instead of the ones compiled in Table IV.
The comparison to phenomenology is difficult, because in our computation fermionic disconnected contributions to the three-point functions have been neglected. However, if one identifies the quantities computed here with an upper limit to what is called valence contribution in phenomenology, we observe that our value for x is smaller compared to phenomenology, while the value for x 2 is also larger compared to phenomenology, but has large error bars.
From the discussion in the previous section it is clear that a computation including fermionic disconnected diagrams is highly desirable. Thus, we are planning to repeat this computation by including fermionic disconnected contributions to the three point functions. Then also the gluonic moments ought to be computed to properly perform the renormalisation procedure. The chiral fit for x gives the following best fit parameters b a -1.6(7) P Z vDD (β = 1.90) 1.31 (2) P Z vDD (β = 1.95) 1.35 (2) P Z vDD (β = 2.10) 1.41 (2) P r (β = 1.90) 5.30 (7) P r (β = 1.95) 5.78 (6) P r (β = 2.10) 7.60 (8) with correlation coefficients in the same order as above: