Valence parton distribution of pion from lattice QCD: Approaching continuum

We present a high-statistics lattice QCD determination of the valence parton distribution function (PDF) of pion, with a mass of $300$ MeV, using two very fine lattice spacings of $a=0.06$ fm and $0.04$ fm. We reconstruct the $x$-dependent PDF, as well as infer the first few even moments of the PDF using the 1-loop perturbative LaMET framework. Our analyses use both RI-MOM and ratio-based schemes to renormalize the equal-time bi-local quark-bilinear matrix elements of pions boosted up to $2.4$ GeV momenta. We use various model-independent and model-dependent analyses to infer the large-$x$ behavior of the valence PDF. We also present technical studies on lattice spacing and higher-twist corrections present in the boosted pion matrix elements.

The above definition involves quark and anti-quark displaced by z − along the light-cone (and made gauge-invariant by the Wilson-line W + (z − , 0) = P exp i z − 0 dz − A + that runs along the light-cone. The dimensionless light-cone distance ν is referred to as the Ioffe-time and the matrix element M(ν, µ), renormalized in the MS scheme by convention, is referred to as the Ioffe-time distribution (ITD). Notwithstanding such a straight-forward definition of PDF, the unequal time separation z − posed a challenge to the Euclidean lattice computation until recently.
To overcome this issue, it was first proposed to use the so-called quasi-PDF (qPDF), which is defined from matrix elements of equal-time bilocal quark bilinear operators and can be related to the PDF for large hadron momenta [3]. This method was then developed into LaMET which provides the framework to calculate all parton physics [4]. Later, there was suggestion to use the so-called pseudo-PDF approach [5,6], which relates * xgao@bnl.gov † nkarthik@bnl.gov the same matrix elements to the light-cone correlations for PDFs at small distances. The hadron matrix element that is central to both LaMET and the pseudo-PDF approaches is E, P z |ψ(z)W z (z, 0)γ t ψ(0)|E, P z ≡ 2E(P z )M(zP z , z 2 , µ R ).
It is very similar to Eq. (1), except that quark and antiquark are at equal-time and separated by spatial distance z and evaluated in an on-shell hadron state at large spatial momentum P z . Thus, zP z can be thought of as the quasi light-front distance measured in units of the Compton wavelength of the hadron and M as the quasi lightfront correlations. Such a matrix element can be easily computed on the lattice [3,4]. In the literature, the matrix element M is also referred to as the pseudo Ioffetime distribution (pITD) [5]. As a crucial step in the UV regulated field theory, the multiplicative renormalizability of the bilocal operator was recently demonstrated to all orders of perturbation theory [7][8][9]. The renormalized matrix element (and its Fourier transforms with respect to z or zP z ) can be systematically related to the PDF within the framework of Large-Momentum Effective Theory (LaMET) [4,10]. The LaMET matching factors from various intermediate renormalization schemes for the equal-time bilocal bilinear matrix element at some renormalization scale µ R to the MS PDF at a factorization scale µ are known to 1-loop accuracy [11][12][13][14][15], and recently, papers related to 2-loop matching have also appeared [16][17][18]. A related good lattice crosssections [19,20] approach has also been recently proposed to calculate PDF on the lattice. In practice, the lattice calculations and the perturbative factors are at fixed order, the different methods may have different advantages and drawbacks. The status of these calculations is summarized in recent review papers [10,[21][22][23].
In this paper, we study the valence pion PDF. The study of pion PDF is interesting for several reasons, both technical as well as with interesting physics issues. The most interesting reason being that the pions are the pseudo-Nambu-Goldstone bosons of QCD and it is important to study its structure in order to understand the relation between hadron mass and hadron structure. Closely related to this, is the question of how fast the PDF vanishes as x approaches 1. This issue of whether the vanishing behavior is (1 − x) 2 or slower is being vigorously debated with various non-perturbative approaches [24][25][26][27][28][29][30][31], now including lattice QCD [32][33][34], thanks to the LaMET formalism. There have been LO and NLO analyses of the experimental data [27,[35][36][37][38][39][40][41][42], but the results are less constrained than the nucleon PDF due to availability of experimental data and therefore, the lattice calculations can have large impact here. The other interesting reasons for studying pion in particular are technical. First, the smallness of the pion mass means that it is easier to have highly boosted hadronic states required in the qPDF approach. Second, the excited state contamination for pions is less problematic due to larger gaps at typical momenta of 1-2 GeV. There has been lattice calculations of pion PDF using the quasi/pseudo-PDF frameworks [34,[43][44][45], and also using the good lattice cross-section approach [32,33].
In our previous work [34], we studied the valence pion PDF in 2+1 flavor QCD using the mixed action with lattice spacing a = 0.06 fm and LaMET approach. In the sea, we used Highly Improved Staggered Quark (HISQ) action, while in the valence quark sector we used clover improved action with hypercubic (HYP) smearing [34]. We extend this study in three ways in this paper. First, we perform calculations at another smaller lattice spacing, namely a = 0.04 fm. Second, we increase the statistics in the a = 0.06 fm ensemble by more than two-fold. Third, we combine the analysis of the bilocal bilinear matrix element renormalized in RI-MOM scheme [46] with the ratio scheme (also referred to as reduced ITD [6]), and also propose and use generalizations of the ratio scheme with the promise of lesser higher-twist contamination. At a practical level, it has been conventional in the lattice calculation that used quasi-PDF formalism to use an intermediate RI-MOM scheme, while those using pseudo-PDF formalism to use an intermediate ratio scheme. We do not make such distinctions, and simply refer to matrix elements of operator in Eq. (2) that is made gaugeinvariant with a straight Wilson-line as bilocal bilinear matrix elements, or simply as the matrix elements for the sake of brevity, in various renormalization schemes; RI-MOM matrix element or ratio matrix element, for example.
The plan of the paper is as follows. In Section II, we discuss the details of the lattice ensembles, statistics and other computational specifics. In Section III, we elaborately describe the extraction of ground-and excited-states of pion from the boosted two-point functions. In Section IV, we describe the extraction of the boosted pion matrix element from three-point function via excited-state extrapolations. In Section V, we dis-cuss the various renormalization schemes used. Readers not interested in the details of the lattice calculation can skip Sections II-V. In Section VI, we describe the twist-2 LaMET formulation which forms the basis of the results presented in the following sections. We also present a study of higher-twist contamination in this section. The Section VII contains the direct extraction of the valence moments of pion from the P z z and z 2 dependences of bilocal bilinear matrix element. In Section VIII, we reconstruct the x-dependent valence PDF at µ = 3.2 GeV based on fits to the pion matrix elements using phenomenology motivated ansatz for the PDFs. In Section IX, we address the issue of large-x exponent of the valence pion PDF based on model dependent fits as well as from a novel model-independent method we introduce here. In Section X, we speculate the continuum results based on our observation at two fine lattice spacings. The conclusion and comparisons with other analyses are given in Section XI. More technical details are present in the appendices.

II. LATTICE SETUP
In this work, we use two different L t × L 3 lattice gauge ensembles both of them with relatively small lattice spacings -(1) ensemble with lattice spacing a = 0.06 fm with lattice extents 48 × 64 3 , and (2) a finer ensemble with a = 0.04 fm with extent 64 × 64 3 . These gauge ensembles were generated by the HotQCD collaboration [47] using 2+1 flavor Highly Improved Staggered Quark (HISQ) action [48] in the sea. In both these ensembles, the sea quark mass was tuned such that the pion mass was 160 MeV. On these gauge field ensembles, we used 1-HYP tadpole improved Wilson-Clover valence quarks. That is, we used the Wilson-Clover quark propagator in the Wick contractions required in the computations of the three-point and two-point functions, and the gauge links that went into the construction of the propagator were smoothened using 1 step of HYP smearing [49]. We set the clover coefficient c sw = u −3/4 0 , where u 0 is the average plaquette with 1-HYP smearing; we used c sw = 1.02868 and 1.0336 for a = 0.06 fm and 0.04 fm respectively. We tuned the Wilson-Clover quark mass m q a in both the ensembles so that the valence pion mass, m π , is 300 MeV. Through an initial set of tuning runs we determined m q a = −0.0388 for a = 0.06 fm and m q a = −0.033 for a = 0.04 fm lattices. For this pion mass, the values of m π L t on the a = 0.06 fm and 0.04 fm lattices are 5.85 and 3.89 respectively. Thus it would be more important to take care of wrap around effects in the finer lattice and we do so in the analysis. With the usage of 1-HYP smeared gauge links in the Wilson-Clover operator, we did not find any exceptional configurations at both the lattice spacings. We used the a = 0.06 fm ensemble in our previous analysis of the valence PDF of pion [34]. With this work, we have increased the statistics used in this ensemble by more than two times. ensemble mqa mπLt nz z range #cfgs (#ex,#sl) a, Lt × L 3 a = 0.06 fm, -0. 0388    The most basic element of this computation is the Wilson-Dirac quark propagator inverted over boost smeared sources and sinks [50] as we discuss more in the next section on two-point functions. We used the multigrid algorithm [51] for the Wilson-Dirac operator inversions to get the quark propagators. These calculations were performed on GPU using the QUDA suite [52][53][54].
We used boosted quark source [50] and sink with Gaussian profile, as we discussed in detail in [34]. Instead of using the gauge-covariant Wuppertal smearing [55] to implement the Gaussian profiled quark sources, we gaugefixed the configurations in the Coulomb gauge to construct the sources as we found it to be computationally less expensive. We fixed the radius of the Gaussian profile on a = 0.06 fm and a = 0.04 fm ensembles to be 0.312 fm and 0.208 fm respectively. We discussed the details of tuning the Gaussian smearing parameters in the Appendix of [34]. Using these quark propagators, we are able to compute hadron two-point and three-point functions in hadrons boosted to momentum P z = 2πn z /(La).
We tabulate the details of the statistics used in the two ensembles in Table I. We increased the statistics in two ways (a) using statistically uncorrelated gauge field configurations, which are labeled as #cfg in Table I, and (b) by using All Mode Averaging (AMA) [56] on each gauge configuration. In order to mitigate the reduction in the signal-to-noise ratio in both the three-point and twopoint functions as one increases P z ∝ n z , we used more gauge field configurations for larger n z than at smaller ones. In a = 0.06 fm ensemble, we effectively increased the statistics 32 times by using 1 exact Dirac operator inversion and 32 sloppy inversions in the AMA per configuration. In the a = 0.04 fm ensemble, we increased the number of exact and sloppy solves for n z = 2, 3 and more for n z = 4, 5. We used a stopping criterion of 10 −10 and 10 −4 for the exact and sloppy inversions respectively.

III. ANALYSIS OF EXCITED STATES IN THE TWO-POINT FUNCTION OF BOOSTED PION
In this section, we discuss the computation of boosted pion correlators and the extraction of the excited state contributions. Using a smeared (s) pion source π s (P, t) for pion π + that is moving with spatial momentum P = (0, 0, P z ) along the z-direction, we computed the twopoint function of pions C ss 2pt (t s ; P z ) = π s (P, t s )π † s (P, 0) .
In this computation, we used momenta on a periodic lattice for n z = 0, 1, 2, 3, 4 and 5 at both lattice spacings. These values of n z correspond to P z up to 2.15 GeV and 2.42 GeV on the a = 0.06 fm and 0.04 fm lattices respectively. For ease of reference, we have tabulated the physical values of P z for the two lattices in Table II. Such large momenta are central to the applicability of the Large Momentum Effective Theory framework. It is important that we are able to suppress the excited state contributions to the two-point function within smaller source-sink separations t s to deal with the signal-to-noise ratio at larger t s . This is the reason for the smeared pion source and sink, π s , that are constructed out of smeared quark fields, u s and d s . We constructed two-kinds of two-point functions: smeared-source (s = S) point-sink (s = P ) correlators referred to as SP, and smeared-source (s = S) smeared-sink (s = S) correlators referred to as SS henceforth. For smeared sources, we used boost smeared Gaussian profiled sources, as is now standard in the lattice PDF computations. We have tabulated the values of the tunable parameter ζ for the boost smearing [50] at different P z in Table II.
The two-point functions enter the PDF determination in two ways; for determining the excited state spectrum of the boosted pion on the two lattices, which in turn will enable us to extract the boosted pion matrix elements. Below, we will discuss the excited state analysis of the two-point function. In our previous publication [34], we discussed the extraction of the pion spectrum in detail for the a = 0.06 fm lattice. Since the only difference in this paper is the increased statistics for this ensemble, we focus on the pion spectrum in the finer a = 0.04 fm lattice in this section. In Fig. 1, we show the effective mass E eff (t s ) of pion at different P z as a function of source-sink separation t s for the SP (open symbols) and SS correlators (filled symbols) respectively. For comparison, the values of E(P z ) for the ground state pion based on its dispersion relation are shown by the horizontal lines. One can notice that the signal-to-noise ratio gets poorer at shorter t s as P z is increased. Therefore, we are forced to work in a range t s /a = 9, 12, 15 and 18 corresponding to physical distances of 0.36 fm to 0.72 fm for the case of three-point functions. The largest operator insertion times τ , which we will discuss in the next section, are t s /2. In this range of t s , the effective mass is not plateaued and careful consideration of excited states becomes important. Up to n z = 3, it is clear that the effective masses plateau at the the dispersion values for the pion. One can also note that SS correlator approaches the plateau faster than SP as expected. The difference between SS and SP correlators is due to the differences in the amplitudes of the states in the two, and we will use this advantageously in the extraction of first and second excited states of the pion.
In order to determine the energy levels E 0 , E 1 , . . ., we fit the spectral decomposition of C 2pt (t s ), A n (e −Ents + e −En(aLt−ts) ), with E n+1 > E n . The above expression is truncated at N state to both the SS and SP two-point function data over a range of values of t s between [t min , aL t /2]. We performed this fitting with one-state (N state = 1), twostate (N state = 2), and three-state (N state = 3) ansatz. As evident from the behavior of effective mass in Fig.  1, in order for the 1-state fits to work, we had to use t min > 0.56 fm and the results were consistent with the one from dispersion relation E 0 (P z ) = m 2 π + P 2 z with m π = 300 MeV. When we performed an unconstrained 4-parameter 2-state fit to both the SS and SP correlators, we found the approach to the expected E 0 (P z ) to be at even shorter t min ∼ 0.2 fm. Since we were able to obtain the ground state energy E 0 (P z ) reliably from one and two exponential fits to both the SS and SP correlators and they agree with the expectation from the dispersion relation well, we then fixed the value of E 0 to its dispersion value to perform a more stable two and three exponential constrained fits with one less free parameter.
The results for the first excited state E 1 (P z ) using different t min in such a constrained two-state fits for n z = 3 and n z = 4 are shown in the top and middle panels of Fig.  2; the top panel is for SP and the middle one for SS. One can notice that for t min /a > 10, it is possible to reliably estimate the first excited state in both SP and SS correlators, and the two estimates are also consistent with each other giving more confidence in the results. The horizontal lines in the figures correspond to the expected result for E 1 (P z ) based on a single particle type dispersion relation E 1 (P z ) = P 2 z + E 2 1 (P z = 0). As t min is increased, the fitted values of E 1 (P z ) are actually the dispersion values. We observed this behavior at different P z as well. We will address this more in the end of this section. Having understood the actual spectral decomposition of the pion correlator, it has been found to be better practically to use the effective value of E 1 and the corresponding amplitude A 1 in the range of τ ≈ t s /2 used in the two-state fits to three-point function [57]. By doing this, we effectively take care of excited states higher than E 1 that could be present at τ ≤ t s /2 in the two-state fits to the thee-point function. We follow this procedure here and take the value of E 1 and A 1 in the pseudo-plateau region for E 1 seen in middle panels of Fig. 2 for t min ∈ [5a, 10a].
We also performed constrained 3-state fits on the SS two-point function. Besides fixing E 0 , we also imposed a prior on E 1 using its best estimate from the SP correlators with the corresponding errors [58]. The results for E 1 and E 2 from this analysis are shown in the bottom panel of Fig. 2 for n z = 3 and 4. As a consistency check, the 3-state prior fit is able to reproduce the input prior for E 1 starting from t min /a = 2. It also results in an estimate for E 2 which is large and noisy, and it is likely that it is an effective third state capturing several higher excited states. For our excited state extrapolations, such an effective estimate is sufficient. We repeated the above set of analysis for the a = 0.06 fm lattice and we were able to obtain the ground and first excited state reliably.
In Fig. 3, we show the first two energy levels for both a = 0.04 fm and 0.06 fm lattices, as a function of P z . It is not very surprising that the ground state, which is the pion, follows the particle dispersion well even up to P z = 2.4 GeV on the fine lattices we use. But, it is remarkable that the first excited state E 1 also follows a single particle dispersion relation. We noted this also in our discussion of Fig. 2. To solidify the claim, we observed the same behavior in both SS and SP channel. Also, the difference between E 1 on the two physical volumes 24.9 fm 3 and 16.78 fm 3 for the a = 0.06 fm and a = 0.04 fm lattices is not seen. Thus, it is likely not a multi-particle state with a gapped finite volume spectrum that mimics a single particle state. In order to account for the 300 MeV pion mass, we added 0.16 GeV to the PDG value [59] of the first pion radial excitation, π 1 (1300) to estimate a value of 1.46 GeV. This value agrees well with our estimates of E 1 (P z = 0) at both the lattice spacings. Therefore, we find it reasonable to conclude that the ground state is the pion and the first excited state is the radial excitation of pion, π 1 , with E 1 (P z = 0) being identified with its mass.

IV. EXTRACTION OF BARE MATRIX ELEMENTS FROM EXCITED STATE EXTRAPOLATIONS
The next ingredient in the extraction of the pion matrix element is the three-point function involving the insertions of smeared pion source π † S (P, 0) and smeared sink π S (P, t s ) separated by an Euclidean time t s and projected to spatial momentum P = (0, 0, P z ). The operator O Γ (z; τ ) is the isospin-triplet operator that involves a quark and anti-quark that are spatially separated by distance z where x = (x, τ ) with τ being the time-slice where the operator is inserted, and the quark-antiquark being displaced along the z-direction by L = (0, 0, 0, z). The operator is made gauge-invariant through the presence of the straight Wilson-line of length z, W z (x + L, x), that connects the lattice sites at x + L to x. The Wilson-line is constructed out of 1-level HYP smeared gauge links to get better signal to noise ratio. The matrix Γ is either the Dirac γ-matrix γ z or γ t for the unpolarized PDFs that we will study in this paper. For the case of lattice Dirac operators that break the chiral symmetry explicitly at finite lattice spacings, it was shown perturbatively in that O γz mixes with the scalar operator O 1 due to renormalization [11,46]. Such mixing is absent in the case of O γt . In addition to this mixing, we also found in our  4. The source-sink ts and operator insertion τ dependence of the ratio R(ts, τ ; z, Pz), at fixed z and Pz = 2πnz/L are shown for the lattice spacing a = 0.06 fm. The top-rows are for nz = 0, middles ones for nz = 3 and bottom ones for nz = 5. The left panels are for z = 0, middle panels for z = 6a and right ones for z = 12a respectively. Each plot has left and right sub-panels. In the left sub-panels, the ts − τ /2 dependence is shown at ts = 8a (red squares), 10a (green circles) and 12a (blue triangles). The corresponding colored bands are the 1-σ errorbands from Fit(2, 3) (see text). On the right sub-panels, the extrapolation (grey band) to ts → ∞ is shown as a function of a/ts at fixed τ = ts/2.
previous work [34] for the case of pion that O γz is comparatively noisier compared to O γt with same statistics, and also suffered from larger excited state contamination. Therefore, we resort to only the usage of Γ = γ t in this paper. The pure multiplicative renormalization of O γt also allows us to explore the renormalization group invariant ratios in addition to RI-MOM scheme as an advantage, and we will explain this in detail in the next section. The above u − d three-point function is purely real in the case of pion, and the real part is symmetric about z = 0. Therefore, we determined the three-point function in both positive and negative z, and averaged over them. In the plots that follow, we will display the three-point function in the positive direction only. In addition, only the quark-line connected piece contributes to the isotriplet three-point function. We refer the reader to the Appendix of [34] for detailed proofs of the above characteristics.
From the three-point function and the two-point function, the central quantity from which the bare matrix element can be obtained from, is the ratio In order to take care of the wrap-around effect due to the finite temporal lattice extent L t , we replace C 2pt (t s ; P z ) with C 2pt (t s ; P z ) − A 0 exp (−E 0 (aL t − t s )) where A 0 and E 0 are the amplitude and energy of the ground state obtained via fits to the two-point function in the last section. This is especially important to take care of at P z = 0 on our lattices. In the above equation, the variables are t s and τ at fixed z and P z , and hence we will keep z and P z implicit in the discussion of R below. Through the spectral decomposition of R, it is easy to see that 1 with E n+1 ≥ E n , E 0 = E π and A n = Ω|π|π . In the infinite t s limit, R(t s , τ ; z, P z ) is equal to the bare matrix element h B (z, P z ) = π|O γt (z)|π . In practice, we obtain h B (z, P z ) by fitting the right-hand side of Eq. The left panels are for z = 0, middle panels for z = 9a and right ones for z = 18a respectively. Each plot has left and right sub-panels. In the left sub-panels, the ts − τ /2 dependence is shown at ts = 9a (red squares), 12a (green circles), 15a (blue triangles) and 18a (pink inverted-triangles). The corresponding colored bands are the 1-σ errorbands from Fit(2, 3) (see text). On the right sub-panels, the extrapolation (grey band) to ts → ∞ is shown as a function of a/ts at fixed τ = ts/2.
the ratio R. The fit parameters are the matrix elements E n , P |O γt (z)|E n , P ) . We take fixed values of E n and A n from our analysis of C 2pt that we discussed in the last section; namely, in two state fits, values of E 1 , A 1 were taken from the pseudo-plateau seen in Fig. 2 that covers the typical range of τ used here, while in the three state fits, the values of E 1 , A 1 were fixed to the actual dispersion values of π 1 and E 2 , A 2 effectively captured the tower of higher excited states. We truncated the number of states N entering the fit ansatz in Eq. (10) at N = 2 and 3. To reduce the excited state contamination, we excluded cases where operator insertion is too close to either the source or sink by using only values of τ ∈ [n sk a, t s − n sk a]. We used n sk = 1, 2 for N = 3 and n sk = 2, 3 for N = 2. We denote such N -state fits as Fit(N, n sk ).
In Fig. 4 and Fig. 5, we show some sample results of the extrapolations using Fit(2, 3) for the a = 0.06 fm and a = 0.04 fm lattices respectively. Each panel in the plot has two sub-panels. Let us first focus on the larger left sub-panels which show the dependence of R(τ, t s ) on τ − t s /2. The lattice data for R are shown as the symbols with the colors distinguishing the different t s . For the a = 0.06 fm lattice, we used t s /a = 8, 10 and 12 (i.e., t s = 0.48 fm, 0.6 fm and 0.72 fm) in the fits. Sim-ilarly, we used t s /a = 9, 12, 15 and 18 for a = 0.04 fm ensemble, which corresponds to similar physical values of t s = 0.36 fm, 0.48 fm, 0.6 fm and 0.72 fm respectively. Along with the data for R(t s , τ ), we have also shown the results from Fit(2, 3) as the similarly colored bands. The result for the matrix element h B , i.e., t s → ∞ limit of the fit, is shown by the grey horizontal band in the figures. The degree to which extrapolation differs from the actual data in the range of t s < 1 fm can be seen from the smaller right sub-panels, where we have shown the 1/t s dependence of the data (points) as well as the fit (grey band) with τ = t s /2, the maximal distance of operator from source and sink. In general, one can see that the extrapolations get steeper as the value of z increases. However, given the small errors at smaller z, the extrapolation again plays a significant role at smaller z. From the agreement of the two-state fits with the actual data, one can gain confidence in the extrapolations.
In addition to the N -state fits, which are sensitive to the values of E n , A n , we also used the summation technique [60] which does not require inputs of the spectral details of the two-point function. For this, we use the For large t s , we would find a linear behavior in t s of R sum as We refer to this method where we ignore O(e −(E1−E0)ts ) corrections and fit only h B (z, P z ) and B 0 as Sum(n sk ). Since our source-sink separations are less than 1 fm, we also included the additional e −(E1−E0)ts correction in the fitting ansatz as We refer to this method as SumExp(n sk ). In Fig. 6, we show a sample result for the summation fits. In the left and right panels of the figure correspond to a = 0.04 fm and 0.06 fm lattice ensembles. We have used momenta P z = 2πn z /(La) with n z = 2 in both the cases at an intermediate separation z = 0.72 fm in both ensembles. The lattice data for R sum are shown as the red circles. The result from a linear fit to the data is shown as the red band. The slope of the fit is the estimator of the matrix element h B . One can see in both the cases that the straight line fit is able to describe the data. However, one can certainly see deviations from the straight line fit at t s = 18a for the a = 0.04 fm case. For comparison, the expectation for R sum (t s ) from the 2-state fit described above is shown as the green band. Here, the curve is able to describe the data at all t s well and can be seen be seen to approach a straight line with larger slope only for t s > 0.72 fm. In order to account for these discrepancies, we also show the result from SumExp as the blue dashed line. This result does deviate from the simple Sum and agrees better with the expected result from Fit. This shows that there are residual O(e −(E1−Eπ)ts ) effects which cannot be ignored in the summation fits in the ranges of t s we are working with. While we have picked an example case where we observe this discrepancy to be larger, similar discrepancy could be seen in other values of P z and z as well in the case of a = 0.04 fm data. The Sum data agreed better with expectation from SumExp and Fit for the a = 0.06 fm data. Therefore, we use the results from Fit, and only use Sum and SumExp to serve as cross-checks on the results.
As we demonstrated above, the t s → ∞ extrapolations lead to values of h B which are not simply obtained from plateau values of R(t s , τ ) even for the largest t s = 0.72 fm we use. Therefore, a way to reasonably justify the correctness of our extrapolations is by adapting the multiple fitting schemes, namely Fit, Sum and SumExp, and show consistency among them. This is what we show in Fig. 7 and Fig. 8, for the a = 0.06 fm and 0.04 fm lattices respectively. The different panels show the results for four different values of P z = 2πn z /(La). In the top part of the different panels, we have shown the bare matrix element h B (z, P z ), obtained by Fit (2,3) as the black open squares, as a function of the length of Wilson-line z. Since we are working with iso-triplet matrix element for the pion, only the real part of h B is non-zero. One should remember that the bare matrix element at any finite lattice spacing has the Wilson-line self-energy divergence, exp(−cz/a), which causes the rapid decay of h B (z, P z ) as a function of z in the figures. With the increased statistics used in our computation, one can note that we are able to obtain matrix elements with good signal to noise ratio even up to momenta corresponding to n z = 5 in both the lattice spacings. Below the top part of each panel in Fig. 7 and Fig. 8, we show the deviations, ∆(z), of different extrapolation methods from values obtained with Fit(2,3) as a function of z. That is, where h B method is the bare matrix element obtained using an extrapolation technique method, which could be Fit(2,2), Fit (3,2), Sum(2), or SumExp(2), in Fig. 7 and Fig. 8. If the extrapolations are perfect, then we would find ∆(z) to be consistent with zero at all z and P z . For comparison, we also show the statistical error in h B Fit(2,3) (z, P z ) as the grey error band along with the values of ∆(z). For a = 0.06 fm case shown in Fig. 7, we find ∆(z) is consistent with zero within error for larger P z while there is little tension at smaller P z in the top two panels. The small but visible deviations of Fit (3,2) is less than 2σ. The deviation of Sum(2) is comparatively larger, but when we supplement Sum(2) with the exponential corrections, i.e., SumExp(2), the ∆(z) moves towards zero and becomes consistent with zero. This again points to the importance of excited state effects that cannot be neglected in summation fits on our lattices. This effect is more apparent in the case of a = 0.04 fm lattice shown in Fig. 8. Thus we understand the deviation of Sum from the rest as an excited state effect, and we find that the Fit(2,3), Fit(2,2), Fit (3,2) and SumExp(2) are all consistent among themselves. Thus, we are able to demonstrate the goodness of our extrapolations. Henceforth, we will use Fit(2,3) for both a = 0.04 fm and 0.06 fm ensembles in discussing our further analysis.
A well determined matrix element that can be used to cross-check our results is the value of bare matrix element at z = 0, which in the continuum limit will be the total isospin of pion, which is 1. At any finite a, the bare matrix element suffers from O(α s (µ = a −1 )) correction to 1, which under finite renormalization will be canceled by Z V . If the excited-state extrapolations were perfect and the finite volume effects were negligible, the estimates of h B (z = 0, P z ) cannot change with P z up to possible finite a corrections at non-zero P z . In order to check for this, we show the behavior of h B (z = 0, P z ) as a function of P z in Fig. 9. For a = 0.06 fm lattice, the value of h B (0, 0) is 1.0404(4) and the values of h B (0, P z ) get smaller than this value gradually at larger P z , albeit only by less than 2% by P z = 2.15 GeV. This P z dependence is likely to arise due to increasing lattice spacing effect at higher momenta, and empirically, it was possible to fit the P z dependence to an ansatz For a = 0.04 fm lattice, the value of h B (0, 0) is 1.045(1) which is higher than value of h B (0, 0) at a = 0.06 fm. However, one expects h B (0, 0) to decrease and approach 1 as a → 0 [61]. One observes a sharp decrease in the value of matrix elements at non-zero P z to values around 1.025 and changes little with P z > 0. We were able to understand this anomalous behavior at P z = 0 to arise from larger periodicity effects (∼ e −Mπ(Lt−ts) ) in the P z = 0 three-point function for the finer a = 0.04 fm lattice (which is in addition to such wrap-around effects in two-point function that we corrected for in Eq. (9)). We discuss this further in Appendix A, and we estimate the value of h B (0, 0) after correcting for the wrap-around effect to be 1.024 (1). For the a = 0.06 fm case, this effect is negligible. The (approximate) corrected estimate for h B (0, 0) is shown as the filled blue square in Fig. 9, which shows surprisingly good agreement with the estimates at other non-zero P z . We used the same fitted (P z a) 2 ansatz that we discussed above, with only the value of a changed from 0.06 fm to 0.04 fm, and the result is shown as the blue dashed curve in Fig. 9. This nice agreement gives credence to our explanation of lattice spacing effect be-  ing the cause of the mild P z dependence in a = 0.06 fm h B (0, P z ) estimates and the even milder P z dependence in a = 0.04 fm estimates. We discuss the estimation of Z V within the RI-MOM framework in Appendix B which give results consistent with the values from the bare pion matrix element in Fig. 9.

V. RENORMALIZATION
The bare matrix element h B (z, P z ) obtained in the last section needs non-perturbative renormalization in order for it to have a well defined continuum limit. The non-perturbative renormalization removes the UV selfenergy divergence of the Wilson-line which is inherently non-perturbative and can only be captured by methods such as the ab-initio lattice QCD heavy-quark potential computations (c.f. Ref [62] for the ensembles used here). With the removal of this non-perturbative piece, one would expect the remaining renormalized matrix element to be describable within the perturbative large momentum effective theory framework. Therefore, a judicious choice of the nonperturbative renormalization scheme for the bilocal quark bilinear operator that is implementable on an Euclidean lattice and at the same time reduces the higher-twist corrections to the matrix element in any given small values of z is important.
RI-MOM is one such renormalization scheme that uses renormalization conditions at off-shell space-like external quark four-momentum P R . A more careful description of the calculation of RI-MOM factor as applied to our work can be found in [34]. The RI-MOM renormalized matrix element is defined as where Z q is the quark wavefunction renormalization factor (c.f. Ref [63]) and Z γtγt is the renormalization factor for O γt (z) defined via the condition imposed using the amputated matrix element evaluated with quark external states at momentum p, Λ(p), as The above condition is referred to as the / p-projection scheme within the RI-MOM scheme [13,46]. The operator O γt does not mix with any other operator, unlike O γz [11,64]. We used the Landau gauge fixed configurations to determine Z γtγt (z) non-perturbatively in both a = 0.06 fm and a = 0.04 fm ensembles. We will refer to the component of P R along the direction of Wilson-line as P R z and the norm of the component perpendicular to zdirection as P R ⊥ . Since the value of h R 0 (z = 0, P z , P R ) = 1  for the pion, we impose this condition through a redefinition This implicitly takes care of the effect of Z q and at the same time reduces the statistical errors in h R at the other non-zero values of z through their correlation with h R (z = 0). Instead of using quark external states, it is possible to cancel the UV divergence in h B (z, P z ) using the pion matrix element at a different fixed reference momentum P 0 z , that is, h B (z, P 0 z ). Such a procedure to remove the UV divergences via renormalization group invariant ratios is referred to as the ratio scheme [10,14]. With this, we can define a renormalized matrix element, The choice P 0 z = 0 has been used in literature and the resulting matrix element M(z, P z , 0) is also referred to as the reduced ITD [6,14]. Non-zero P 0 z was applied to proton in [57]. Similar to the RI-MOM matrix element, we can reduce the statistical errors by redefining the matrix element as so that the condition M(z, P z , P 0 z ) = 1 is automatically fulfilled. We use values of P z > P 0 z in this work. The preference for using P z , P 0 z > Λ QCD will become clearer  10. Comparison of renormalized matrix elements at fixed Pz = 1.45 GeV in the ratio scheme with generalized non-zero values of the reference momentum P 0 z . The different colored symbols correspond to different P 0 z = 0, 0.48 and 0.97 GeV. The effect of changing P 0 z is significant, but it does not cause a big difference in signal-to-noise ratio at smaller z that we are interested in. It is the expectation that matrix elements with Pz, P 0 z > {ΛQCD, mπ} suffer from lesser highertwist contamination.
with the discussion on perturbative matching in the next section. In Fig. 10, we compare the result of M(z, P z , P 0 z ) for three different P 0 z = 2πn 0 z /(La) for n 0 z = 0, 1 and 2 on a = 0.06 fm lattice. These values of n 0 z correspond to 0, 0.48 and 0.97 GeV respectively, and thus using even the lowest n 0 z available makes sure P 0 z > Λ QCD . The effect of using P 0 z as a new scale leads to significant changes to the P z z and z 2 dependence, which will be taken care of the corresponding twist-2 expressions. But one should note that we do not significantly compromise on the quality of signal by choosing non-zero values of P 0 z < 1 GeV, and hence, they are as good choices of the the reference momentum scale in the ratio scheme as P 0 z = 0 GeV. Our choice of the normalization conditions in Eq. (17) and Eq. (19) such that the value of pion matrix element at z = 0 is 1, assumes implicitly that our estimates of the matrix elements at z = 0 do not suffer from any systematic corrections. In the discussion around Fig. 9, we found about 1% systematic errors at z = 0 due to deviations of the matrix element as a function of P z . Below, we justify that the imposition of the normalization conditions Eq. (17) and Eq. (19) also reduces some of these systematic errors. Instead of imposing the normalization multiplicatively as in Eq. (19), an equally good choice is additively through The multiplicative and additive normalization are equivalent, only provided M 0 (0, P z , P 0 z ) is itself exactly 1. In Fig. 11, we compare the result of M add (z, P z , P 0 z ) and M(z, P z , P 0 z ) at P z = 1.29 GeV on a = 0.04 fm lattice. The left and right panels are for P 0 z = 0 and 0.48 GeV respectively. For comparison, we have also shown the matrix element M 0 before imposing the normalization. First, one can note the error reduction due to the normalization at smaller values of z. As we discussed in the last section, the z = 0 matrix element at P z = 0 for a = 0.04 fm suffers from larger systematic effects than the rest. From the left panel which shows the result for P 0 z = 0, we surprisingly find that the difference between M add and M is absent within the errors at all z. On the right panel, which uses P 0 z = 0.48 GeV, the agreement is perfect between all the estimates of M. Through this, we demonstrated that the systematic effects in our matrix element determination are further reduced due to the ratios using the prior knowledge that the local matrix element at z = 0 is 1 for pion.
Finally, we address the lattice corrections to the renormalized matrix elements. In Fig. 12, we have shown the comparison of renormalized matrix elements at two lattice spacings a = 0.06 fm (red circles) and 0.04 fm (blue squares) plotted as a function of ν = zP z . The top and bottom panels show the comparison using RI-MOM and n 0 z = 1 ratio scheme respectively. Due to lattice periodicity constraints, we could only choose pion momentum P z that are approximately the same at the two lattice spacings; namely, P z = 1.29 GeV for a = 0.06 fm and P z = 1.45 GeV for a = 0.04 lattices. By looking at the pion matrix elements at the two lattice spacing as a function of P z z, such small mismatch between P z should affect results only logarithmically in this discussion. For the RI-MOM scheme, we have chosen a comparable set of renormalization momenta (P R z , P R ⊥ ) = (1.93, 2.23) GeV for a = 0.06 fm lattice and (1.93, 2.51) GeV for a = 0.04 fm lattice. In the bottom panel, we have used matrix element in ratio scheme with n 0 z = 1. We find only a little difference between the matrix elements at the two fine lattice spacings. To aid the eye, we have also shown bands that cover ±2% variation on the a = 0.06 fm and a = 0.04 fm data. Within this band, the real parts of the data are consistent with perhaps little more correction to the imaginary part of RI-MOM at intermediate ν. Thus, we can bound the lattice corrections in our data to be at the level of 1 to 2%. In the RI-MOM data, perhaps there are residual lattice spacing effects of about 1% at different z. Even though the this lattice spacing effect is only about a percent, we will see that a 2 P 2 z corrections become important in the analysis at smaller z due to to their very small errors ensured by the normalization process. The computation of renormalized pion matrix element is the final step as far as the non-perturbative lattice input is concerned. The perturbative matching lets us make the connection between the renormalized boosted hadron matrix element with the light-cone MS PDF, f (x, µ). Since the renormalization factors for the RI-MOM h R (z, P z , P R ) and the ratio M(z, P z , P 0 z = 0) do not depend on the PDF of the hadron itself, they lead to simpler factorized expressions and hence let us consider them first. Using such expressions, we will consider M(z, P z , P 0 z ) for non-zero P 0 z . Taking Ji's proposal [3,4] of quasi-PDF in the RI-MOM scheme,q(x, P z , P R ), which is the Fourier transform of the z-dependent matrix element the perturbative matching is expressed as a convolutioñ The kernel of the convolution C RI (x, . . .) is of the form (with dependence other than x being implicit), where C (1) (x) is the 1-loop contribution [11][12][13][14], and the notation [. . .] + represents the standard plus-function 2 .
Though we have used RI-MOM scheme in the above equations, one can use the matrix element in ratio scheme, M(z, P z , P 0 z = 0), as well with a corresponding C ratio (x, P z , µ).
An equivalent approach, that is suitable for our analysis in the real space z, instead of performing a Fourier transform in Eq. (21) to the conjugate x, is through the formulation of operator product expansion of the renormalized boosted hadron matrix element [14] using only the twist-2 operators. That is, for the case of M(z, P z , P 0 z = 0) computed at large P z and with z in the perturbative regime, its OPE that is dominated by twist-2 terms is Here, x n (µ) are the nth moments of the PDF at a factorization scale µ, The coefficients c n (z 2 µ 2 ) are the perturbatively computable Wilson-coefficients defined as the ratio of MS Wilson-coefficients, c n ( The 1-loop expressions for c n (z 2 µ 2 ) can be found in Refs. [14]. These Wilson-coefficients are related to the matching kernel C ratio through the relation [14], The corrections denoted as O z 2 Λ 2 QCD arise from the operators in the OPE that are of twist higher than two.
For the RI-MOM scheme, a similar OPE that is valid where the RI-MOM Wilson-coefficients are c RI n . Using the multiplicative renormalizability of the bilocal operator O γt , we can deduce that (28) where Z ratio→RI is the perturbatively computable P zindependent conversion factor from RI-MOM to ratio scheme [11,21]. By taking the ratio of Eq. (26) for the ratio scheme and a corresponding similar expression for the RI-MOM scheme involving c RI n and C RI , we can work out the conversion factor Z ratio→RI to be up to 1-loop order. Though there is an explicit P z present in the above expression, its dependence gets canceled in the final expression, as expected.
We can now consider the ratio scheme for general values of P 0 z . Noting that M(z, P z , P 0 z ) = M(z, P z , 0)/M(z, P 0 z , 0), we can write the twist-2 expression as Such an expression cannot be written in a factorized form involving a convolution of a perturbative kernel and PDF. As we noted in the beginning of this section, we anticipated this since the "renormalization factor" is (M(z, P 0 z , 0)) −1 for the ratio scheme at non-zero P 0 z and hence by itself dependent on the hadron PDF, unlike the RI-MOM or the P 0 z = 0 ratio schemes. However, as far as the practical implementation of the analysis is concerned, the non-factorizability of Eq. (30) is not a hindrance, and the analysis proceeds in exactly the same way for all the schemes considered, i.e., by extracting the moments x n from the boosted hadron matrix elements either in a model independent way or by modeling the PDFs to phenomenology inspired ansatz. The reader can refer to Ref [57] for this method implemented for the nucleon.
Finally, the above discussion ignored any presence of lattice spacing corrections present at smaller z at the order of few lattice spacings that could spoil the applicability of the twist-2 expression as it is. As discussed in Section IV, we found indications of (P z a) 2 corrections to matrix element at z = 0. Such lattice corrections were removed at z = 0 by taking the ratio and making z = 0 renormalized matrix elements to be one by construction. However, such a procedure will not ensure cancellation of (P z a) 2 corrections at any non-zero z. We will take care of such correction by including fit terms, ra 2 P 2 z , by hand in the twist-2 expressions above, with r being an extra free parameter. As a concrete example, we will modify Eq. (30) to to accommodate for any short-distance lattice artifacts. It is easy to see that the effect of such a (P z a) 2 correction is to shift the second moment in (z/a) −2 manner, in all the twist-2 expressions above. Indeed, we will present evidence for the presence of such (P z a) 2 corrections, and we defer that discussion to Section VII. One should note that the above ansatz for correcting (aP z ) 2 effects is strictly true only for z > 0, since the ratio has to be exactly 1 at z = 0. The actual form of lattice correction would automatically ensure this, but we found this simpler form to be practically enough to describe the z > 0 data, starting from z = a.
B. Numerical investigation of higher-twist effect in pion matrix element at low momenta In the remaining part of this section on perturbative matching, we discuss a way to use the hadron matrix elements at smaller momenta to understand the importance of higher twist effects at intermediate values of z ∼ 0.3 fm to 1 fm, and thereby, understand the rationale for the ratio scheme which hitherto had been discussed using a conjectured factorization of higher twist effects [6]. The Eq. (24), Eq. (27) and Eq. (30) are valid only up to O(z 2 Λ 2 QCD ) higher twist effects. At any large value of P z , such O(z 2 Λ 2 QCD ) are much smaller compared to any of the twist-2 terms. This is the basis of the large momentum effective theory. As a corollary, the matrix element where the higher twist effects show up significantly is the P z = 0 matrix element. This is not a useful observation when applied to the ratio M(z, P z , P 0 z = 0), for which c 0 has the value 1 at all z when P z = 0. This agrees with the twist-2 expectation by construction, but the corrections 1.08 The phase (top) and the magnitude (bottom) of the Pz = 0 matrix element in RI-MOM scheme with (P R z , P R ⊥ ) = (1.29, 2.98) GeV are shown. The expectations from the 1-loop leading twist results are shown as the red bands. In the bottom panel, the actual lattice data is shown using filled black circles. The absolute part clearly suffers from a leading z 2 correction shown as the blue straight line. The twist-2 target mass correction is shown as the green curve for comparison. The lattice data after subtracting the z 2 correction term is shown using open circles.
could show up at other non-zero P z . Therefore, we use h R (z, P z = 0, P R ) in the RI-MOM scheme for this study where we compare the lattice result with the non-trivial z dependence from twist-2 term. Also, since the wraparound effects in P z = 0 matrix elements are negligible only for the a = 0.06 fm lattice, we use this case for this study.
For P z = 0, the only non-zero twist-2 contribution is from the local current operator because all other terms have explicit factors of P z and they become zero. Its z dependence comes from the Wilson-coefficient c RI 0 (z), which is the conversion factor Z ratio→RI (z, P R , µ). We use 1loop expressions to calculate Z ratio→RI . We vary the scale of α S (µ) that enters Z ratio→RI from µ/2 to 2µ with µ = 3.2 GeV, and gives an estimate of the expected error on perturbative result. It is convenient to separate Z ratio→RI and the lattice result h R (z, P z = 0, P R ) into their magnitudes and phases. The phase Arg h R (z, P z = 0, P R ) is the same as Arg Z γtγt (z, P R ) , which is a property of the RI-MOM scheme itself. On the other hand, the magnitude |h R | depends on the pion matrix element.
In the top panel of Fig. 13, we compare the z 2 depen-dence of the phase Arg h R (z, P z = 0, P R ) with the perturbative twist-2 phase Arg Z ratio→RI (z, P R ) . We have chosen a renormalization scale (P R z , P R ⊥ ) = (1.29, 2.98) GeV on the a = 0.06 fm ensemble as a sample case, but the observations hold for other cases as well. We find a good agreement within the perturbative uncertainties up to 0.7 fm, and the lattice data slightly overshoots the 1-loop result for larger z. Nevertheless, the overall qualitative agreement validates the 1-loop perturbation theory as applied to quark external states used in RI-MOM Zfactor. This should serve as a companion observation to the studies on RI-MOM Z-factor presented in our previous work [34].
In the bottom panel of Fig. 13, we compare the z 2 dependence of the magnitude h R (z, P z = 0, P R ) with Z ratio→RI (z, P R ) . The actual lattice data is shown as the filled circles. It is clear that the non-perturbative result disagrees with the near constant behavior of the twist-2 term at larger z, and that this disagreement comes from a striking z 2 dependence at larger z. The coefficient, k, of the z 2 dependence is −(63 MeV) 2 (with little variations around this value with P R ), and thus it is reasonable to identify such a term to arise from a higher twist operator or an effective contribution of a number of higher twist operators. There could also be corrections to the leading twist result coming from the twist-2 target mass correction (TMC) [65,66] (and discussed in Appendix C). The 1-loop result with TMC is shown as the green dashed line in the bottom panel of Fig. 13, which is visibly small compared to the observed discrepancy. Numerically, the coefficient of z 2 from twist-2 target mass correction term is −m 2 π x 2 v /8 = −(35.2 MeV) 2 , which about one-third of the observed value (assuming x 2 v ≈ 0.11 as we will see later). In addition, when we correct for the z 2 effect by subtracting it from the lattice data, shown by the open black circles in the figure, we find a nice agreement with the 1-loop, twist-2 expectation. It is quite remarkable that such a simple O(Λ 2 QCD z 2 ) effect is enough to describe the nonperturbative data even up to 1 fm. Now, we take the hypothesis that the observed z 2 effect in h R (z, P z = 0, P R ) is the dominant higher twist effect, and try to understand its effect on the matrix element in the ratio scheme. Perturbatively, the ratio scheme is defined via the subtraction of the UV divergence by a division with n = 0 MS Wilson coefficient, c MS 0 . On the lattice, we identify this procedure as the division by P z = 0 matrix element, and hence the equality in Eq. (24). The underlying assumption is that the higher twist effect in P z = 0 matrix element is negligible or somehow cancels with the higher twist effect present also in the non-zero P z matrix elements. In order to understand this, we can redefine the ratio scheme that better agrees with the assumptions that go into LaMET; namely form the ratio after subtracting off the higher-twist effects where k is the coefficient we determined using the analysis of h R (z, P z = 0, P R ) and assume the same kz 2 correction is present at non-zero P z as well. For k = 0, M = M. The result of this improved ratio M (z, P z , 0) is compared with the usual ratio M(z, P z , 0) in the left panel of Fig. 14 for the first non-zero momentum P z = 0.43 GeV on a = 0.06 fm lattice. The difference between the two ways of defining the ratio are consistent within errors, with perhaps very little difference at larger z. This provides a better understanding of how the O(Λ 2 QCD z 2 ) corrections in the numerator and denominator of Eq. (33) almost cancel each other without resorting to any factorization of higher-twist corrections, and instead, results simply from the smallness of k. Having demonstrated the inconsequential role of higher twist effects in M for z < 1 fm given the errors in the data, we now look closely at the RI-MOM h R (z, P z , P R ) at the same small momentum P z = 0.43 GeV. Within the twist-2 framework, we can obtain h R from M via h R (z, P z , P R ) ≡ Z ratio→RI (z, P R )M(z, P z , 0). (34) In the right panel of Fig. 14, we compare h R , shown as the red band, with h R which are the black filled symbols. We find a deviation from the twist-2 expectation h R for z > 0.3 fm. When we correct for the z 2 effect using |h R (z, P z , P R )| − kz 2 , shown as the open symbols, we find a very good agreement with h R . Putting together the above results, we self-consistently justified that the observed kz 2 effect in P z = 0.43 GeV is almost the same as in P z = 0 as we assumed, and that M is least affected by such corrections. At higher momenta P z , such highertwist effect will play even lesser role for z < 1 fm. The blue curve is the expectation for the zPz dependence at z = a based on the moments obtained at z = 5a without correcting for (Pza) 2 lattice terms. The red curve is obtained after accounting for such lattice corrections.

VII. A MODEL-INDEPENDENT COMPUTATION OF THE EVEN MOMENTS OF VALENCE PION PDF
In this section, we apply the LaMET formalism, that we discussed in Section VI, to our lattice data for the isovector u − d PDF of pion. For this, we will use the boosted pion matrix element in the ratio scheme with non-zero reference momentum P 0 z = 2πn 0 z /(La) with n 0 z = 1 and 2. This way, we expect to suffer from smaller non-perturbative corrections and also avoid the larger periodicity effect in zero momentum matrix elements (we also discuss results using n 0 z = 0 for a = 0.06 fm lattice, where wrap-around effect was small, in Appendix H). Through Eq. (30), we can find the values of the moments x n by fitting them as free parameters such as to best describe the zP z and z 2 dependence of the M(z, P z , P 0 z ) data. Such a method, usually referred to as OPE without OPE, has been previously applied to the case of reduced ITD (n 0 z = 0) for pion [44] and nucleon [67][68][69].
A. Connection between isovector PDF and valence PDF of pion First, it is important to recall as to how the u − d isovector pion matrix element that we compute on the lattice relates to the valence PDF of pion. Let f u (x) and f d (x) are the u and d quark PDF with support in x ∈ [−1, 1] and a convention that includes the quark distribution for x > 0 and anti-quark distribution for Due to isospin symmetry in The moments that occur in the OPE expressions Eq.
This could be understood from the fact that u parton could only be produced radiatively in π + , and hence f u only has a sea quark distribution, which thereby cancels the sea quark distribution of f u in the above definition. Due to the isospin symmetry present in our QCD computation that does not include QED corrections, . Unlike the u − d PDF of pion, both even and odd valence moments x n v of the pion are non-vanishing. By comparing the above equivalent definition of the valence PDF in terms of u and d quark PDFs in pion with the u − d PDF in Eq. (35), one can deduce that and that for the moments x n v , n is even, 0, n is odd.
Thus, the OPE expression in Eq. (30) for u − d pion matrix element has only even powers n, from which we could obtain the values of x n u−d for even values of n, which as we discussed is the valence moment x n v . Unfortunately, the u − d matrix element does not directly let us access the odd valence moments, but we will later try to determine them based on models of valence PDF f π v itself.

B. Method for model independent fits
We performed model independent determinations of x n v by fitting the rational functional form in Eq. (30), which we denote as M OPE (z, P z , P 0 z ) here, with the even moments x n v as the fit parameters, over a range of z 1 ≤ z ≤ z 2 and P 0 z ≤ P z ≤ P max z . The possibility of larger lattice corrections at very short separation , z 1 , has to be accounted for. Therefore, we tried fits including or excluding (P z a) 2 correction term in M OPE (z, P z , P 0 z ) as discussed in Section VI. For a larger range [z 1 , z 2 ], there is a larger curvature in the data for M, which makes the fits sensitive to the higher order terms of ν in Eq. (30). On the other hand, by using a larger z 2 , there is the undesired possibility of working in a nonperturbative regime of QCD. We strike a balance between the two by choosing the maximum, z 2 , over range of values from 0.36 fm to 0.72 fm. We choose the factorization scale µ to be 3.2 GeV in the following determinations. Since the Wilson coefficients c n are known only to 1-loop order, the scale of strong coupling constant α s is still unspecified. We take care of this perturbative uncertainty by using the variation in Eq. (30) when the scale of α s is changed from µ/2 to 2µ as part of error, where µ is the factorization scale at which x n v are determined. Concretely, we minimize the following χ 2 to determine the moments: While the above expression is a convenient way to include the perturbative error in the analysis, it comes at the cost of missing the covariance matrix. We take care of it by using the same set of bootstrap samples for all z and P z . We use the factorization scale µ = 3.2 GeV to determine α s used in the twist-2 expressions; for this, we used the values α s = 0.33, 0.24 and 0.19 at scales µ/2, µ and 2µ respectively, by interpolating the running coupling data compiled by the PDG [59]. Since we take the variation of α s with scale into account in the error budget of our analysis, a precise input of α s is not necessary. We can also improve the estimate of higher moments by imposing priors on N prior lower moments by using where x i prior and σ prior i are the prior on i-th moment and error on the prior respectively. We used this method only to determine x 6 v with prior imposed on only x 2 , or both x 2 and x 4 v . For the prior, we used the result of fits with z 2 = 0.5 fm and the error on that estimate as σ prior . In the future, it would be interesting to use estimates of lower moments from the other twist-2 localoperator techniques on the same gauge ensemble as priors in the LaMET methods in order to determine higher moments.
We point out an improved way to implement the fit for valence pion PDF. Naively, one might expect that including more terms in Eq. (30) will lead to unstable results and larger errors due to the increase in the number of fit parameters, x n v . For the case of valence PDF pion, we can use an additional fact to constrain the moments -that of the positivity of f π v (x), and hence of This stems from the fact that the u-quark is present at the order O(α 0 s ) due to its valence nature while the d-quark is only in the sea and hence its distribution can start only at O(α s ). Thus, it is a well justified expectation that The positivity of f π v (x) leads to the conditions that the even derivatives of x n v with respect to n are positive (i.e., d m x t v dt m | t=n = x n log m (x) v > 0 for even m) and that the odd derivatives (i.e., m is odd) are negative. The interesting consequences are the inequalities These inequalities lead to strong constraints on the fitted moments and lead to the stabilization of the estimates (and their errors) of the lower moments as one increases the number of terms in Eq. (30) to larger values, thereby eliminating the order of Eq. (30) as a tunable parameter and prevents over-fitting the data. The two inequalities in Eq. (41) can be easily implemented through a change of variables where the sum runs over even i and j for the pion. The parameters λ j > 0, and N being the largest even moment used in the fit. In the discussions below, we used even moments up to x 8 v in the fits over multiple z. In cases where certain higher moments were irrelevant to the fits, they promptly converged to values very close to zero without affecting the relevant smaller moments. In this way, we do not have to choose the order of the polynomial to be used in the fits.

C. Determining an estimate, its statistical and systematic error
Since the various estimates in this section and the rest depend on the range [z 1 , z 2 ] and the value of n 0 z used, we define the central estimate of a quantity A and its systematic error as Mean(A) and SD(A) respectively; here, Mean(A) is the mean over different estimates (variations in fit range etc.,) in a given bootstrap sample, and SD(A) is the standard deviation of various estimates of A within the same bootstrap sample. The notation Mean(A) and SD(A) stand for average of those mean and standard deviation over the bootstrap sample. In this way, we obtain the statistical error on Mean(A) also in the standard bootstrap procedure. We will use this procedure in the later sections too, and the extra dependences on model ansatz, and renormalization schemes (ratio, RI-MOM scheme and their various scales) will also enter in evaluating the systematic error.

D. Model independent analysis of moments at fixed z 2
The ν = zP z dependence can come from either the z variation at fixed P z or from P z variation at fixed z. We first look at the latter case. In Fig. 15, we show the result of fitting the rational polynomial in ν given by Eq. (30) to the a = 0.06 fm data at different fixed z. For the case shown, P 0 z = 0.43 GeV. In this analysis at fixed z, we did not take any (P z a) 2 correction into account. The results of the fits at various fixed z are shown as the bands having the same color as the corresponding data points. Since we only have five different values of P z , the smaller z data cover shorter ranges in ν compared to the larger fixed z data. It is clear from the data that in order to be sensitive to deviations from simple ν 2 term, we need to resort to data at larger z > 8a = 0.48 fm as well. We repeated this analysis with n 0 z = 2 also. In Fig. 16, we show the value of x 2 v that is extracted from the fits as a function of the fixed values of z used in the fits. The results as obtained from both n 0 z = 1 and 2 are shown in the left and right panels. In order to look for lattice spacing effects, we have shown results from the two lattice spacings (but keeping in mind that the two n 0 z at two lattice spacings lead to slightly different P 0 z in physical units). The inferred moments are more precise for n 0 z = 1 than for n 0 z = 2, as one would expect from deteriorating signal as momentum is increased. One can see a plateau in x 2 v starting from 3a even up to z = 0.7 fm. This shows that the z 2 dependence in the pion matrix element is canceled to a good accuracy by the perturbative Wilson coefficients c n (µ 2 z 2 ).
There is a clear tendency for the fitted x 2 v to increase at very short lattice distance z/a ∼ 1 which is most likely a result of increased lattice corrections at smaller z. One can see this by comparing the results from the two lattice spacings and noting that at fixed short physical distance z, there is tendency for the a = 0.04 fm data to lie closer to the plateau than the a = 0.06 fm data. If the lattice spacing effect is coming from (P z a) 2 corrections, then we should find the (z/a) −2 behavior of x 2 as we outlined in Section VI (where we ignore the logarithmic dependence in z present in c 2 for this first analysis in this subsection.) The fits to such x 2 + 2r(z/a) −2 are shown as the corresponding colored curves in Fig. 16. Indeed, we find a very nice description of the observed data, thereby, show the importance of (P z a) 2 corrections at first few lattice separations z/a. Also, as a consistency check, the values of r from the fits on the two lattice spacings were about the same, namely, 0.021 and 0.022 on the 0.04 fm and 0.06 fm lattice spacings respectively. It could be counter-intuitive to find a rather large lattice spacing effect affecting the moments when we do not find anything unusual about the small z/a in Fig. 15, or in Fig. 12 where we compared the data at two different lattice spacings. In order to understand this, we take the values of moments as obtained from z = 5a (which lies in the plateau of x 2 v ) and reconstruct the expected ν dependence at fixed z = a using Eq. (30), without including any (P z a) 2 corrections in the expression. In the inset of Fig. 15, we compare this expected curve (blue) with the actual z = a data points. The clear disagreement between the two is the cause of the anomalously large x 2 v ≈ 0.15 at z = a in Fig. 16. One should note the rather enlarged scale on the y-axis of the inset, and the disagreement is actually subpercent. But, the data at small z/a is so precise that such small lattice spacing effects show rather clearly in the extracted moments. This is the crux of the problem. After accounting for the ra 2 P 2 z correction, the expected curve is shown in red, which agrees perfectly with the data and gives x 2 v that is consistent with the one extracted from larger z/a. In the analyses henceforth, we will use the correction term r(P z a) 2 term in the fits as outlined in Section VI with r being an extra fit parameter, and this way, we were able to use z 1 = a, 2a, 3a in the fits and obtain no contradictory strongly z 1 -dependent results.

E. Model independent combined analysis of moments
In order to estimate x n v , it is better to to fit both zP z and z 2 dependence using all the data within z ∈ [z 1 , z 2 ] and P z > P 0 z . In Fig. 17, we show the best fit values of x 2 v as a function of the maximum of the range of z, i.e., z 2 . The left and the right panels are for the a = 0.06 fm and a = 0.04 fm data. Along with z 2 dependence, we have also shown x 2 v from the three different values of z-range minimum, z 1 = 1a, 2a and 3a (as we noted, we include a r(P z a) 2 term in the fits in order to be able to use z 1 = a and 2a). The two different colored symbols differentiate the reference momenta n 0 z = 1 and 2. These combined fits with moments being the fit parameters lead to typical χ 2 /dof ≈ 0.7 in all the cases. For both the lattice spacings, we find the various estimates to be consistent with each other. The scatter of values at a = 0.04 fm seems to be centered around a slightly lower value than at a = 0.06 fm, pointing to a small lattice spacing dependence. Using the convention for summarizing the various estimates, we find at µ = 3.2 GeV, with the first error being statistical and the second one being systematic. We input the fit results from z 1 = 1a, 2a, 3a, z 2 ∈ [0.24, 0.6] fm, and n 0 z = 1, 2 to obtain the above single estimate. These estimates with statistical error band, and with both statistical and systematic error band are shown in Fig. 17. For comparison, the estimate of x 2 v from JAM collaboration [70] at µ = 3.2 GeV is at a slightly lower value, 0.095. The softgluon resummed ASV result [27] is even lower at about 0.086 at the same scale µ.
In Fig. 18, we show a similar plot for x 4 v at µ = 3.2 GeV. At each z 2 , we show determinations with z 1 = 1a, 2a and 3a. We find consistent determinations with various fit ranges and renormalization procedures. We estimate These estimates are the bands in Fig. 18. The JAM estimate, x 4 v = 0.032, is slightly lower than for the 300 MeV pion studied here [70]. Whereas the ASV result for the fourth moment can be inferred to be about 0.023. For both x 2 v and x 4 v , we did not use priors. On the other hand, it was not possible to obtain a good estimate of x 6 v without inputting the knowledge of the lower moments using the procedure we outlined previously. We obtained results for x 6 v at the two lattice spacings by inputting prior for only x 2 v , and by using priors for both x 2 v and x 4 v . We display the results for the latter case in Fig. 19. In addition, we have to use z 2 > 0.5 fm in order for x 6 v to be a relevant parameter in the fit. As we noted in the previous section, the Λ 2 QCD z 2 corrections seem to be canceled effectively even in the ratio scheme with P 0 z = 0, and the error we commit by using values of z 2 up to 1 fm might not be large and also further reduced by non-zero P 0 z we use in the modified ratio scheme. Perhaps this is the reason, we find the estimates to be independent of z 2 and P 0 z to a good degree. We estimate To compare, the JAM and ASV estimates as inferred from their fits are 0.015 and 0.009 respectively. In the above fits, we obtained the coefficient r of the (P z a) 2 correction to be −0.026(7)(10) and −0.018(8)(8) for 0.06 fm and 0.04 fm lattice spacings, which are quite consistent with each other as expected, and with our rough estimate in the last subsection. We should also point out that in the above discussion, we did not include any target mass correction (trace) terms in the OPE used in fits since we did not find any significant change by including such additional terms due to the smallness of pion mass.

VIII. VALENCE PDF OF PION BY FITS TO BOOSTED PION MATRIX ELEMENTS IN REAL SPACE
In the last section, we estimated the even moments directly from the equal-time boosted pion matrix elements. However, it is not possible to reconstruct an x-dependent PDF using only the knowledge of the first few even moments. One way of PDF reconstruction from the boosted pion matrix element is through data interpolation over the range of z where lattice data is available and then extrapolate it to zero smoothly at larger z [15,71]. Instead, as in our previous work, we adopt the method of using phenomenology motivated ansatz for f π v (x) and fit the ansatz to our lattice matrix element over ranges of z smaller than 1 fm. In this way, we avoid the usage of data with z 1 fm which could be deep in the nonperturbative regime, and might not be consistent with the perturbative framework that we rely on. There are also other methods of PDF reconstruction that have been investigated in the literature [72,73].

A. PDF ansatz and analysis method
As is typical in the global analysis of valence PDFs, we use two different valence pion PDF ansatz (46) with the first one being a special case of the second and hence more restrictive. The normalization factors N , N are chosen such that 1 0 f π v (x)dx = 1. The parameters α, β, s, t are the tunable fit parameters. These model PDFs enter the analysis via their corresponding moments, for example x n (α, β) = 1 0 x n f π v (x; α, β)dx, which appear in the OPE expressions; Eq. (30) for the ratio scheme and Eq. (27) for RI-MOM scheme. In both the schemes, we corrected for (P z a) 2 lattice artifacts that affect smaller z by using a term r(P z a) 2 in the OPE expressions with r being a fit parameter, as we did in our model independent fits. Through this, we can construct the model matrix elements M model (z, P z , P 0 z ; α, β, . . .) and h R model (z, P z , P R ; α, β, . . .). Let us first consider the ratio scheme. In addition to the statistical error σ stat (z, P z , P 0 z ) for the lattice data point M(z, P z , P 0 z ), there is also the perturbative uncertainty resulting from the 1-loop truncation of the twist-2 Wilson coefficients. We quantify this error through the arbitrary nature of the scale µ of the strong coupling α s (µ), as we did in Section VII. We use µ in the α s to be the same as the factorization scale of the PDF, and quantify the error we commit through the systematic error σ sys (z, P z , P 0 z ) which we define as the change in M model (z, P z , P 0 z ; α, β, . . .) when α s is changed from α s (µ/2) to α s (2µ). That is, Let us take the JAM data and the ASV analysis data at µ = 3.2 GeV as a specific case. The JAM data can be described to a very good accuracy by the form Eq. (46) with α = −0.37 and β = 1.20. In Fig. 20, we show the result for M model (z, P z , P 0 z ) at P z = 1.29 GeV, P 0 z = 0.43 GeV using the JAM valence PDF [70] with solid curves, and using ASV result [27] using dashed curves. For each case, we plot three different curves for M model as obtained using α s (µ/2), α s (µ) and α s (2µ). For comparison, the actual lattice data and the error for M model (z, P z , P 0 z ) is also shown. We can see that the spread in M model for both JAM and ASV get especially important for z > 0.4 fm, and become comparable to the statistical error in the data. Therefore, given the significant perturbative uncertainty that is unavoidable at present, it would be misleading to favor or rule out models of PDF (such as JAM and ASV results in the example here) simply based on the statistical precision of the lattice data. Therefore, we select the model PDFs that best describes the shape of the lattice matrix element that takes σ sys into account, by minimizing, The correlations between the lattice data at different z and P z are partly taken into account by picking M(z, P z , P 0 z ) from the same bootstrap samples. Similarly, in the case of RI-MOM matrix element, we fit only the real part, Re h R (z, P z , P R ) , and the imaginary part is obtained as an outcome. In the case of RI-MOM matrix element, we found taking care of σ sys to be even more important as the α s dependence starts from c RI 0 , unlike in ratio scheme.
B. Results for f π v (x) In Fig. 21, we show the resulting best fit model matrix elements along with the actual lattice data. We have used the 4-parameter ansatz f π v (x; α, β, s, t) for the fits shown.
For the fits, we used µ = 3.2 GeV in the perturbative Wilson coefficients, and hence the model PDF corresponds to this factorization scale. In the results shown in Fig. 21, we used only the boosted matrix elements with the quarkantiquark separations z ∈ [2a, 0.5 fm], but we performed the analysis also with z 1 = a, 2a, 3a and z 2 ∈ [0.36, 0.72] fm. The matrix elements at different fixed P z are differentiated (by their color and symbols). In the top and bottom panels we have shown the results for a = 0.06 and 0.04 fm respectively. We have shown the results for the ratio scheme with P 0 z = 2πn 0 z /(La) for n 0 z = 1 and 2 in the left and middle panels of Fig. 21. For n 0 z = 1 ratio scheme, we used the momenta with n z = 2, 3, 4, 5, and for n 0 z = 2, we used n z = 3, 4, 5. For the RI-MOM scheme, we used only the larger set of momenta corresponding n z = 3, 4, 5. This is to avoid the larger O(Λ 2 QCD z 2 ) corrections in the RI-MOM scheme observed in Section V. The results of the fit to the RI-MOM matrix elements at renormalization scales (P R z , P x µ = 3.2 GeV; a = 0.06 fm was used for this reconstruction of valence PDF via their ability to describe pion matrix elements in real space in different ratio schemes involving ranges of quark-antiquark separation z ∈ [2a, 0.5] fm. In each panel, such results for f π v (x) based on ratio schemes with reference momenta P 0 z = 2πn 0 z /L with n 0 z = 1, 2 are shown as different colored bands. For comparison, the JAM result at the same µ is shown as the black band.
the cases gave good χ 2 /dof between 0.5 and 1, and we discuss this in Appendix D. We refer the reader to Appendix H for a similar discussion on fits to P 0 z = 0 ratio matrix elements (i.e., reduced ITD).
Each of the best fit model matrix elements in Fig. 21 correspond to valence PDFs, f π v (x; α, β, s, t) at µ = 3.2 GeV. In Fig. 22, we have shown the results of the valence PDFs, f π v (x; α, β, s, t), that are reconstructed from M(z, P z , P 0 z ). The left and the right panels of Fig. 22 show f π v (x) and xf π v (x) as functions of x. The red and blue bands are for the two values of n 0 z = 1, 2 respectively. For comparison, the JAM valence PDF [70] at the same µ is shown as the black band. At a qualitative level, it is reassuring that the PDFs we determined compares well with the phenomenological result. At both lattice spacings, the results from different P 0 z differ only by a little, and such variations belong in the systematic error budget. However, when we look closely, one can find that the best fit PDFs always have a tendency to be above the JAM result for x > 0.6. This ties back to the PDF moment determination in the last section where we found x 2 v and other higher moments also to be consistently higher than the phenomenological result. In Fig. 23, we show similar results for PDF as obtained using the RI-MOM h R . The results using two different renormalization scales P R are consistent with each other as one would expect. One can also note that the RI-MOM results also agree overall with the one from ratio scheme. When we focus on specific details of the PDF, as we would do next, the difference across renormalization schemes and renormalization scales will become easier to notice.     24. Comparison of PDF as extracted using a simple two-parameter x α (1 − x) β ansatz and from a four-parameter For the results we discussed above, we limited ourselves to a specific fit range in z from 2a up to 0.5 fm using an ansatz f π v (x; α, β). The obvious addendum to this discussion is to also specify what happens when we change the various choices we used in the fits. First, we check that the constructed PDF is not sensitive to the PDF ansatz. We used both the ansatz in Eq. (46) in our analysis, and in fact, the simpler ansatz f π v (x; α, β) by itself is sufficient to describe our pion matrix elements in real space; in all the cases χ 2 /dof varied between 0.5 to 0.9. The ansatz f π v (x; α, β, s, t) includes terms that affect only the small-x behavior and therefore more flex-ible. In Fig. 24, we compare the best fit PDFs using the two ansatz for a sample case that used ratio scheme with n 0 z = 1. It is clear that the ansatz dependence is very little, and the effect of including more free parameters in f π v (x; α, β, s, t) is to increase the uncertainties in the fitted PDFs without changing the overall shape.
It would be cumbersome to describe one dependence after another in terms of the resulting PDFs. Therefore, we summarize the results of fitted PDFs using various choices for z-range [z 1 , z 2 ] and the renormalization schemes via their first four moments Fig. 25. It is noteworthy that even though we cannot access the odd moments directly, we can obtain them indirectly from model PDFs. Let us focus on one of the panels in Fig. 25 to unpack the details. Each point is an estimate of the moment labeled in the x-axis of the panel. The red and blue points result from using f π v (x; α, β) and f π v (x; α, β, s, t) respectively, and demonstrates the variation due to fitted ansatz. For each of the ansatz (i.e., red or blue), the variation due to the range of z used, [z 1 , z 2 ], is shown as one moves up along the y-axis. The results with the same [z 1 , z 2 ] are enclosed within the dashed lines. For each [z 1 , z 2 ], the variation coming from the renormalization scheme used for the equal time matrix elements is shown. We have shown four such renormalized results with each set of [z 1 , z 2 ] -ratio scheme with n 0 z = 1, 2 (denoted as ratio-1 and ratio-2 in the figure), and RI-MOM scheme at two different (P R z , P R ⊥ ) (denoted as RI-1 and RI-2). For a = 0.06 fm, the two RI-MOM scales are (1.29,2.98) GeV and (1.93,2.98) GeV, and for a = 0.04 fm, the two scales are (1.93,3.34) GeV and (2.9,3.34) GeV. It is satisfactory that the results for x 2 v and x 4 v obtained here indirectly agrees well with the direct determination in the last section, and serves as a cross-check. One sees comparatively larger renormalization dependence in the more stringent two-parameter ansatz, but it seems to be reduced and accounted for by using a more flexible four-parameter ansatz. As we change z 2 from 0.42 fm to 0.72 fm, the results remain almost intact. From these fits, we estimate the moments and their statistical and systematic errors (coming from [z 1 , z 2 ], RI-MOM and ratio renormalization schemes, and the two PDF ansatz) as The indirect determination of the first moment x v shows that each of the two valence quarks carry about a quarter of the pion energy as has been seen before. Especially here, one certainly sees a lattice spacing effect that tends to make x v closer to the JAM value of 0.223 at µ = 3.2 GeV. Such a lattice spacing effect is seen to a lesser extent in x 2 v , and difficult to see in the higher moments.

IX. DISCUSSION ON LARGE-x BEHAVIOR
The large-x behavior of PDFs are of the form characterizing how f (x) vanishes in the x → 1 limit. The exponent β is hadron dependent, and one uses the Brodsky-Farrar quark counting rule [74] to find the typical value of β for a hadron; β = 2 for the pion and 3 for the unpolarized nucleon valence PDFs respectively from this counting rule. In fact, for the proton, one does find the value of β to be close to 3 for the u-quark valence PDF, and motivates the usage of quark counting rule to predict the value of β for other hadrons. But, the value of β from the analysis [37] of E-0615 Fermilab data was found to be about 1 as a contradiction to the quark counting rule. The importance of soft gluon resummation in the analysis of DIS data close to x → 1 limit was pointed in [27], and consequent reanalysis of Fermilab result suggested a value of β ≈ 2. The recent global Monte Carlo analysis of experimental data from JAM collaboration [70] suggests β = 1.2, and concurred by another analysis using xFitter [75]. Nevertheless, quark counting rules are not direct predictions of nonperturbative QCD and there have been lot of recent works on computing β for the pion that relies on alternative nonperturbative arguments, such as DSE, BSE and light-front quantization methods; many such recent attempts [24,25] suggest a value β ≈ 2, but some [28][29][30]76] of them suggest values close to 1, or a cross-over from β=1 to 2 behavior very close to x = 1 [26]. Thus, the issue of the value of β for f π v is still not settled and the lattice computations, as the present one, can play an important role.

A. Model dependent estimate of β
The recent lattice computations [32][33][34] of pion PDF, including our previous work, have attempted to address the issue of β based on the assumption of ansatz of the type in Eq. (46) for f π v (x). In a similar way, we summarize our results on the large-x exponent β and the small-x exponent α in Fig. 26 based on the model dependent analysis that we presented in the last section. The notation and the arrangement of data points in Fig. 26 is the same as outlined in Fig. 25 for the moments. From the plots for β in both the lattice spacings, we find that the fits prefer a value of around 1, and sometimes even smaller than 1. As suggested in [33], the usage of the 4-parameter ansatz does lead to somewhat larger values of β than obtained using the 2-parameter ansatz, but these values are still closer to 1. Even though the moments that correspond to these fits showed little renormalization scheme dependence, the exponents themselves show a larger sensitivity to the lattice renormalization scheme used, with a tendency for the RI-MOM scheme to consistently give lesser values of β compared to those from ratio schemes. One can also notice a somewhat increasing tendency of β when larger fit range [z 1 , z 2 ] is used. It is possible that the favored value of β could be slightly larger than our estimates if we were to include data at larger values of zP z , but we have restricted z 2 to be less than 0.72 fm to remain close to the perturbative regime. A naive argument would suggest the LaMET formalism does not permit one to access smaller values of x < Λ QCD /P z , we find that the model dependent analysis clearly gives robust values of α ∼ −0.5. This is because the pion matrix elements constrain the few low moments, which in turn are functions of both α and β due to the model used. To summarize, the overall shape of the PDF and its first few moments are well determined by the usage of phenomenology motivated ansatz, but the exponents α and β themselves show sensitivity to the renormalization schemes as well as the range of z used. Nevertheless, the estimates of large-x exponents from this analysis have a tendency to lie closer to 1 rather than 2. Quantitatively, we estimate taking into account the different fit ranges and ansatz dependences in the ratio schemes as well as RI-MOM scheme. These are the bands shown in Fig. 26. Since FIG. 25. The first four valence PDF moments x n as inferred from the best estimates of PDFs f π v (x) that best describes the equal-time pion matrix elements in ratio and RI-MOM renormalization schemes. The top panels are for a = 0.06 fm and bottom ones for a = 0.04 fm. The dependence of x n on all the variables in the lattice analysis is summarized in the above plot. The foremost variable is the range of quark-antiquark separations used in the fits z ∈ [z1, z2]. Such variations are bunched together as blocks separated by the dashed lines along the y-axis. The second variable factor is the renormalization scheme of the matrix elements: it could be RI-MOM scheme or ratio scheme at reference scale P 0 z . At fixed z ∈ [z1, z2], four different renormalization points are shown: ratio scheme at n 0 z = 1 (ratio-1), n 0 z = 2 (ratio-2), RI-MOM scheme at two different scales P R , denoted as RI-1 and RI-2. This scheme and scale variations are bunched together within the dotted lines. The tertiary variable is the fit ansatz: the results obtained using the ansatz f π v (x) = N x α (1 − x) β are shown in red and those using RI-MOM scheme has a tendency to obtain smaller β systematically, and since we found the ratio scheme performs better at suppressing the higher-twist effects, we also give the estimates below using only the n 0 z = 1 and n 0 z = 2 ratio schemes Indeed, leaving out RI-MOM scheme results leads to slightly larger β, but still around 1. For comparison, the JAM global fits at the same µ give α = −0.37 and β = 1.20.
The downside of the above model dependent analysis is the question of whether by using a sufficiently general functional form f π v (x), it is possible to find β ≈ 2. For example, if we performed the analysis with β = 2 fixed, the χ 2 /dof was between 1.5 and 2 as opposed to the global minimum between 0.5 and 1 when β was allowed as a free parameter. We discuss such an analysis at fixed β = 2 in Appendix I. A recent lattice study [33] using the good lattice cross-section approach found β = 2 is not ruled out when an ansatz, which we refer to as the 4-parameter ansatz here, is used while β ≈ 1 was preferred when the 2-parameter ansatz was used. Another recent study [77] . .) that best describes the ratio and RI-MOM real space data. The first two panels are for a = 0.06 fm and last two for a = 0.04 fm. The dependence of the exponents on all the variables in the lattice analysis is summarized in the above plot. The notation is similar to Fig. 25. The plot shows the effective large-x exponent β eff (n) as a function of n. The data points are obtained by using values of moments x n v obtained from the model independent combined fits. The smaller error bar is only the statistical error and larger error bar is statistical plus systematic error. The red and blue points are the results using our a = 0.06 fm and 0.04 fm respectively. The red and blue bands are the expectation for the behavior of β eff (n) as obtained from the model-dependent fits to the pion matrix elements. The green curve is the expectation using the JAM data [70] and the purple dashed curve is obtained using the ASV analysis [27].
found that with the limited sensitivity of the lattice calculations to higher moments, it is difficult to make definite conclusions about the large-x behavior. Therefore, we discuss a novel model independent way to find β.

B. Model independent estimate of β
We note that the higher moments get more contribution from larger x, and hence, are more sensitive to the exponent β. Consequently, one finds that the moments x n approach zero in the large-n limit in a manner dependent only on β as The exponent is universal, and independent of the smallx (i.e., x α ) or intermediate-x (i.e., G(x)) behaviors, but the constant of proportionality in Eq. (53) does depend on the details of the PDF. We outline a proof of this behavior in Appendix F. The asymptotic behavior of largen moments was also considered in the context of evolution of β with scale in Refs [78,79]. Thus, one can determine β in a model independent way by taking the log-derivative of the above behavior, A discretised form of the above expression that is suitable for a practical implementation is by defining an effective value of β at finite n as As one uses moments at larger values of n in the above equation, one will find β eff (n) to plateau at the value of large-x exponent β. While this method is straight forward, it also points to the challenge of addressing the large-x exponent -one needs to compute larger moments for a reliable estimate, and puts a limit on what lattice studies can actually address about β without a modeling bias.
In Section VII, we determined the first few even moments in a model independent manner. We used these estimates of x 2 v , x 4 v , x 6 v and x 8 v from the model independent analysis using Eq. (55). Since, the larger moments are required we used the values from the fits where prior was imposed on x 2 v . Even though the larger moments are relatively noisier, the ratios of moments that enter Eq. (55) are better determined owing to their correlations. We estimated the central values of β eff and its statistical and systematic error by the outlined procedure to take care of the variations in [z 1 , z 2 ] and P 0 z . We show the result of β eff (n) as a function of n in Fig. 27 from the two lattice spacings (red and blue data points). We notice that it is possible at the most to use data up to β eff (n = 6). As one would expect, the higher order 1/n corrections to the n −β−1 behavior to be the largest for n = 2, and hence, we find β eff (n = 2) to have larger value around 3.5. For n = 4 and n = 6, the value of β eff decreases and stays around 1.5; the errors are large enough to see any n dependence beyond n = 2. Thus taking the well determined estimate at n = 4 as a proxy for the β, we find Thus, we find some evidence for β to be between 1 and 2, consistent with both our model dependent findings of β ≈ 1 and with quark-counting rule expectation of 2. Thus, we are unable to rule out β = 1 or 2 simply from this model independent analysis. Also, apriori, it is not clear what large-n means; whether one will observe an approximate plateau for β eff (n) at n ∼ O(1) or O(100). Using our model dependent fits as well as JAM result, we will show some evidence below that the plateau is likely to develop for n ∼ O(1) for pion.
In order to see the expected behavior of β eff for n > 6, we simply use x n (α, β, s, t) from our PDF fits in the last section to compute the corresponding β eff (n). In this case, we know that β eff (α, β, s, t) → β in the large-n limit, and we already found such model dependent analysis predict β ≈ 1. We show this resulting β eff (n; α, β, s, t) as the red and blue bands in Fig. 27. We used the parameters as obtained from combined fits to M(z, P z , n 0 z = 1) using a fit range [z 1 , z 2 ] = [a, 0.72 fm] for the case shown. We also show the expected result for β eff using the JAM result as the green curve. The important observation here is that the models of the type in Eq. (46) predict that β eff (n) is almost plateaued by n = 4, and makes β eff (n ≥ 4) to be meaningful estimators of β. To contrast with the JAM expectation for β eff , we also plot β eff as expected using ASV soft-gluon resummed analysis as the purple dashed curve (in order to infer the higher moments for the ASV result, we interpolated their result evolved to µ = 3.2 GeV with The β eff for ASV never goes below 2, and approaches its plateau value at 2.48 from below. C. A semi-model-independent analysis of pion matrix element and exponent β Based on the asymptotic behavior of large-n moments, we propose a new way to fit the moments to zP z and z 2 dependence of pion matrix elements and, at the same time, obtain the value of β in a manner that is not dependent on PDF ansatz. We fit low moments up to an order N asym in the usual manner and use the asymptotic expression for the moments beyond the order N asym using Eq. (53) with some 1/n corrections, as The fit parameters are the lower moments a 2 , a 4 , . . . , a Nasym−2 , and the parameters β, A 0 , A 1 , A 2 that model the large-n moments. We input the constraint that a 2 > a 4 > . . . > a Nasym−2 > x Nasym (β, A 0 , A 1 , A 2 ). Using this model for the moments in Eq. (30), we fit the parameters to best describe M(z, P z , P 0 z ) in the same way as we described in Section VII, but we use z up to 0.72 fm in this analysis. Some analysis bias comes from the choices of N asym and the order of 1/n corrections to use. We used N asym = 2, 4 and 6 in our fits, and any usage of more than A 2 in the fits made the fits unstable and did not converge properly (it is an asymptotic series after all). It was quite surprising that we were able to even use N asym = 2 in the analysis to get good fits, i.e., using the hypothesis that moments starting from x 2 v can be described by the asymptotic expression in Eq. (57). The resulting values of x n v using this method compares well within errors with the moments obtained in Section VII and Section VIII. For example , using the ratio scheme with n 0 z = 1 in the a = 0.06 fm lattice and with N asym = 4, we get x 2 v = 0.110(3), x 4 v = 0.038(5), x 6 v = 0.016(4) which compares well with the other determinations (and for this case, the other parameters were β = 1.20 (21), A 0 = 0.87 (21), A 1 = −0.51 (47), A 2 = −0.28 (28) and χ 2 /dof = 46.5/48). The novel outcome of this analysis is the estimate of β in addition to moments, and for the n 0 z = 1 ratio matrix element using fits from z = a to 0.72 fm, we get , for a = 0.04 fm. (58) As one relaxes the order N asym where large-n asymptotic behavior sets in from N asym = 4 to 6, the best fit values of β changes from a smaller value ≈ 1.0(2) to ≈ 1.8 (5). Thus, this analysis suggests that the values of β around 1 seem preferred when one is aggressive on the order N asym , but β is consistent within 1-σ error (albeit a noisier estimate) if one uses conservatively larger N asym ≥ 6. We reach the same conclusion again; in order to obtain a conclusive result on β, we need even more precise data at larger P z z to be sensitive to higher moments.

X. CONTINUUM ESTIMATES
In the last part of the paper, we discuss the continuum estimates of the PDF and its moments. The estimates are speculative because we only have two lattice spacings, nevertheless both very fine. We should note that we already demonstrated the presence of lattice spacing effects of the type (P z a) 2 at distances of the order of few lattice spacings and took care of them in our analysis. Once the (P z a) 2 artifacts were removed, it was possible to describe the boosted pion matrix elements at any finite lattice spacing using the twist-2 OPE expressions. Therefore, we assume that any additional lattice spacing effects will simply affect the values of the extracted moments themselves. That is, we model the moments x n (a) at any fixed lattice spacing a to behave as where x n v is the continuum value and d n are numerical coefficients that can be fit to the data. We repeated all the analysis (model-independent estimates of even moments, fits to model PDFs, the semi-model dependent analyses) presented in the previous sections using combined fits to both a = 0.06 fm and 0.04 fm data (for fixed physical values of z 2 , and keeping n 0 z to be the same between the two lattice spacings) using the above ansatz for x n v (a) in the twist-2 OPE expressions; concretely, we used fits of the type M(z, P z , P 0 z ; a) = n c n (z 2 µ 2 ) x n + a 2 d n (−iPzz) n n! + r(aP z ) 2 n c n (z 2 µ 2 ) ( x n + a 2 d n ) for the ratio schemes with P 0 z = 0 for this analysis. We found using d 2 and d 4 as additional fit parameters was sufficient to describe the data at both the lattice spacings with χ 2 /dof between 0.5 and 1 depending on the range of fits and analysis type. Since we found that the ratio scheme succeeded in reducing higher-twist effect well, we focus on this scheme in order to discuss our best estimates and their continuum estimates. As a sample result from this analysis, in Fig. 28, we show the a 2 extrapolation of x 2 v and x 4 v as obtained from the above analysis at fixed z 2 = 0.5 fm and z 1 = a using n 0 z = 1 for both lattice spacings. For the case shown, we used the 4-parameter PDF ansatz for the combined fit. The two data points in the plot are the values for the same case at fixed a = 0.04 fm and 0.06 fm. We remind the reader that this is not a straight-line fit to the two data point, but rather an outcome of the combined analysis as described above (with d 2 = 3.1(1.8) fm −2 , d 4 = 0.79 (62) fm −2 , and χ 2 /dof = 60.5/92.) In Table III, we tabulate all our estimates from different kinds of analysis from the previous sections using only the ratio schemes with n 0 z = 1, 2, where we expect the higher-twist corrections to be even milder. Along with the estimates at two fixed a, we also tabulate our continuum estimates based on the above analysis for each quantity. There is only little effect from including a 2 corrections. We find that removing the lattice spacing effect have a slight tendency to bring the moments closer to the JAM values of 0.223, 0.095, 0.052, 0.032 for the first four moments. The best fit values of the large-x exponent β from the fits, however, continue to remain closer to 1, and thus, the lattice spacing effects might not be an issue in our results.

XI. CONCLUSIONS
In this paper, we presented a lattice computation of the MS isovector u−d parton distribution function of 300 MeV pion and its moments using the recently proposed Large Momentum Effective Theory (LaMET) framework. Using isospin symmetry, we related the properties of isovector pion PDF to the valence u − u PDF of pion, π + .
In order to access the short distance physics required for the LaMET framework, we used two lattices ensembles with very fine lattice spacings of a = 0.06 fm and 0.04 fm for the first time in such pion PDF computations. Using high statistics, we were able to compute the required equal-time bilocal quark bilinear matrix elements evaluated with pions boosted up to 2.42 GeV. Thus, a major advancement resulting from this work is the demonstration that current lattice calculations can satisfy both the theoretical requirement of sub-fermi separations (in order to be consistent with the OPE-based framework reliant on naive power counting for operator hierarchy), and the requirement of large hadron momentum (in order for a controlled truncation of the OPE TABLE III. Summary of results from analyses presented this paper at µ = 3.2 GeV. Each row is a type of analysis, namely -(a) the model-independent estimates of the moments (b,c) fits to 2-and 4-parameter PDF ansatz, (d) a semi-model independent analysis based on modeling x n v by the asymptotic formula for n ≥ 4, and (e) the estimate of exponent β from β eff (n = 4). The columns are the outcomes; namely, the value of the first four valence moments x n v , the parameters α, β, s, t in the PDF . For each analysis, values from two different lattice spacings a are given, and also our continuum expectations, denoted by a → 0, based on a 2 extrapolation are also given. For these estimates, we used the ratio scheme with n 0 z = 1, 2.
at leading-twist that underlies the LaMET). As a handle on quantifying perturbative uncertainties and other higher-twist systematics, we used multiple renormalization schemes for the equal-time matrix element, namely RI-MOM, ratio scheme and new variants thereof with the advantage of reducing higher-twist effects. As a technical elaboration, we proposed and used the pion matrix element at zero pion momentum as a suitable quantity to study higher-twist effects and demonstrated practically as to why the ratio renormalization scheme effectively eliminates higher-twist effects to a good accuracy (with respect to typical errors in lattice data at larger z) even up to 1 fm distances.
From the renormalized boosted pion matrix elements, we performed two kinds of analysis. In the first kind, we obtained the first few even valence moments x n v by fitting both z 2 and P z z dependence of the matrix elements making use of perturbative matching coefficients at 1-loop order. Though the LaMET methodology helps us access higher moments without the problem of mixing, especially with the usage of priors on lower moments, it comes at a cost of introducing dependencies on the range of z and P z z used in the fits, and we discussed them in this work. Folding in such dependencies in the systematic error, we estimated the MS moments x 2 v , x 4 v and x 6 v at µ = 3.2 GeV. In the second kind of analysis, we reconstructed the x-dependent valence pion PDF by modeling the PDF via x α (1−x) β G(x) type ansatz, and fitting the parameters of the model so as to best describe both z and zP z dependence of pion matrix elements in various renormalization schemes with z restricted to subfermi values. We summarize our reconstructed valence PDFs at µ = 3.2 GeV from the two lattice spacing in the top panel of Fig. 29 -the estimate using only statistical error is shown as darker band, and statistical and systematic error (coming from fit ranges, renormalization scheme used, and the PDF ansatz used) is shown as the lighter band. This model dependent method lets us access the odd moments of valence PDF as well. We also provided estimates of the values of moments as well as the PDF in the continuum limit based on the results at the two lattice spacings; the numerical results from various analysis approaches are summarized in Table III. In the bottom panel of Fig. 29, we have shown our estimate for the PDF in the continuum limit for the 300 MeV pion.
We discussed the large-x behavior using our modeldependent PDF that we reconstructed from fits. We found that even though the overall x dependence of the reconstructed PDFs remained the same with variations coming from the fit range, renormalization scheme and PDF model used, the specific details such as the value of the large-x exponent β showed a larger dependence on such analysis choices, but largely showed a tendency to be close to 1. To avoid such issues, we proposed a new model independent observable, constructed out of moments, that converges to the large-x exponent β as one uses larger moments. At present, our computed matrix elements are sensitive only to moments up to order 4 or 6, and given this limitation, we find the effective value of the exponent to be between 1 and 2, but ruling out neither. However, it is at present hard to conclude if such an estimate would remain unchanged as one includes larger moments, which would be possible with increased statistics in the future. But, with this work, the computation that one should perform to reach a model independent robust conclusion about β is clear.  Finally, we compare our PDF determinations with other global fit analysis in the summary plot in Fig. 29. Along with our determinations of valence PDFs at the two lattice spacing (top panel) and our estimate of valence PDF in the continuum limit (bottom panel), we have also shown the Fermilab E-0615 estimate [37] (green symbols), the ASV reanalysis of the Fermilab result after taking soft-gluon resummation into account [27] (dashed green line), the recent JAM Monte-Carlo global analy-sis [70] (black band) and the result from analysis using xFitter [75] (purple line), all evolved to the same scale µ = 3.2 GeV. One can find an overall agreement of our determinations with the phenomenological results; with better agreement with JAM, xFitter and the initial E-0615 estimates, than with the ASV result. Some caveats are clear -firstly, our computation is for 300 MeV pion, and hence a future computation with physical pion mass is crucial. Secondly, we used 1-loop matching coefficients to match the lattice results to MS PDF, and it is at present unclear what the effect of adding higherloop perturbative terms in the matching kernel (and also ASV-type resummation of soft-gluon contribution in the matching kernel, if at all possible) on the extracted PDFs and moments will be (very recently, works [16][17][18] related to 2-loop matching appeared as the present manuscript was being completed). We leave these questions for the future. The value of m π L t measures quantitatively the magnitude of wrap-around effects due to lattice periodicity in the temporal direction. Its value on the 48 3 × 64 lattice is 5.85 and on 64 4 lattice is 3.9. Thus, we expect the wrap-around effect in correlation functions to be more in our finer lattice than in the coarser lattice. We take care of this wrap-around effect in two-point function by replacing A 0 exp (−m π t) with A 0 [exp(−m π t) + exp(−m π (L t − t))]. But there are additional effects in the three-point function itself. To quantify the effects of finite L t on the three point function we recall that finite L t can be interpreted as inverse temperature, and therefore, we can write with H being the QCD Hamiltonian. Inserting the complete set of energy eigenstates we can also write the above expression as C 3pt (t s , τ ) = m,n,k m|π|n n|O γt |k k π † m e −τ En e −(ts−τ )E k e −(Lt−ts)Em . (A2) We can split the above sum into two parts; namely, a part without any wrap-around effect in which the state |m is the vacuum |0 , and the states |n and |k run over excited states with quantum number of pion. The second part captures the wrap-around effect, and in which case, the state |m is the pion state, while the states |n and |k run over states with vacuum quantum numbers. In the discussion below, we restrict these states to include |0 and the first lightest iso-singlet, G-parity positive state, |S , with energy E S . That is, we write the spectral decomposition of the three-point function as C 3pt (t s , τ ) = n,k∈iso-triplet 0|π|n n|O γt |k k π † 0 e −τ En e −(ts−τ )E k + π|π|0 0|O γt |S S π † π e −(ts−τ )E S e −(Lt−ts)Eπ + π|π|S S|O γt |S S π † π e −tsE S e −(Lt−ts)Eπ .
The terms in the sum within the brackets is the part without wrap-around effect, and we used this part in the main text to extract the matrix element. The two terms below that are due to finite L t . We focus on O γt (z = 0) now, where we can make wellmotivated estimates of the wrap-around effect. In this case, the term in the second line in Eq. (A3) involving the vacuum vanishes. Further, we make the following assumptions in order to estimate the wrap-around effect: 1. The state |S is a two-pion state |π, π with both the pions with zero relative momentum, and projected to be iso-singlet, G-parity positive.
With these assumptions and using Eq. (A3), the ratio R = C 3pt /C 2pt becomes R(z = 0, t s , τ ) ≈ π|O γt (z = 0)|π 1 + 2e −EπLt 1 + e −EπLt +(t s , τ )dependent excited state terms. (A4) Thus, we have to correct our estimated value for π|O γt (z = 0)|π from excited state fits by the factor above. For non-zero P z , the values of E π are large and the correction factor is almost 1. For zero P z , this wrap around effect is the highest. For a = 0.06 fm lattice, this factor is 1.0028, which almost unity. However, for the finer a = 0.04 fm lattice, the correction factor at zero P z is 1.020, which is comparable to the estimated value of matrix element at z = 0 itself, and hence, it cannot be neglected. In Section IV, we estimated π|O γt (z = 0)|π = 1.045(1) without taking wrap-around effect in the three-point function into account. We estimate the corrected value to be π(P z = 0)|O γt (z = 0)|π(P z = 0) = 1.024(1), (A5) for a = 0.04 fm lattice. This is comparable to other values of z = 0 bare matrix elements at non-zero P z for a = 0.04 fm lattice (refer Fig. 9). Thus, we understand quantitatively the underlying issue in P z = 0 matrix element for the finer lattice, and hence we avoided the usage of it in the analyses discussed in the main text.

Appendix B: Discussion on ZV
In the main text, we normalized the z = 0 renormalized pion matrix elements to 1, thereby avoiding the issue of vector current renormalization factor. Here, we provide details of the renormalization constant Z V of the vector current O µ =ψγ µ ψ in the RI-MOM scheme for our two lattices. We can write Z V = Z γtγt Z q , where Z γtγt is the renormalization of the vertex function for γ µ = γ t and Z q is the renormalization of the quark field. We use the same notation as in Ref. [34], where the details of RI-MOM renormalization can be found. In Fig. 30, we show Z V for the two lattice spacings used in our study, a = 0.04 fm and a = 0.06 fm as function of the renormalization point p R . To minimize the discretization effects, the lattice momenta are substituted by p µ = sin(ap µ ), so p 2 R = µ=1,4 (p µ ) 2 . The vector current renormalization constant should not depend on p R , because in the a → 0 limit the local current is conserved. Nevertheless, we see a significant dependence on p R . This dependence can be caused by lattice artifacts as well as by non-perturbative effects that for large values of p R can be parametrized by local condensates. As we use off-shell quark states in Landau gauge in the RI-MOM renormalization procedure the lowest dimension local condensate is the dimension two gluon condensate A 2 [80,81]. Lattice artifacts show up as breaking of the rotational symmetry on the lattice. We see from Fig. 30 the fish-bone structure in the lattice data at the level much larger than the statistical errors on Z V . All these effects need to be taken into account if one wants to extract Z V . These effects are easier to understand by analyzing Z q and Z γtγt , separately as discussed below.
In Fig. 31 and 32, we show the numerical results Z q and Z γtγt as function of p R . The numerical results on Z q have much smaller statistical errors compared to Z γtγt . The relative statistical errors on Z q are always smaller than 4.5 · 10 −4 and for large p R are smaller than 5 · 10 −5 . We parametrize the p R dependence of Z q and Z γtγt by the following form Z i = Z 0 i + B/(ap R ) 2 + C · (ap R ) k 1 + C 4 ∆ (4) + C 6 ∆ (6) , i = q, γ t , where This form is motivated by the 1-loop lattice perturbation theory [82,83] and the perturbative analysis with dimension two gluon condensate [84]. For the non-perturbative clover action k = 2, while for Wilson action k = 1. For HISQ smeared clover action with tapdole improved value of c sw we expect O(a) discretization errors to be proportional to α 2 s with a very small coefficient, so it is reasonable to assume that the dominant cutoff effects scale like a 2 . Nevertheless, we also perform fits using k = 1. In Ref. [83] the condensate contribution was ignored but it was included in the analysis of PNDME collaboration [61] when fitting the p R dependence of the renormalization constant. For Z γtγt and a = 0.06 fm the fits with k = 1 and k = 2 work well, while for a = 0.04 fm the fits with k = 2 work better. Fits of Z γtγt with k = 1 give χ 2 /df that is around 2 for a = 0.04fm . It is obvious from Figs. 31 and 32 that the condensate contribution to Z γtγt is quite small, while it is large for Z q .
The fits of Z q have very large χ 2 /dof, most likely because of the very small statistical errors. The value of the condensate obtained from the fit is compatible with the value g 2 A 2 ∼ 4 GeV 2 found in Ref. [81] for a = 0.06 fm. For the smaller lattice spacing it is, however, is twice larger, which could be due to instabilities in the fits. Therefore we fix the condensate to the above value in order to stabilize the fits. From the fits, we obtain the values of Z 0 i which can serve as estimates for Z q and Z γtγt for the two lattice spacings. Multiplying these two renormalization constants we obtain Z V . The results of our analysis are summarized in Table IV. The large uncertainties in Z q come from the differences in the fits with the condensate contribution being fixed and treated as the fit parameter. We can compare our result for Z V at a = 0.06 fm with the value Z V = 0.945(15) from the PNDME collaboration. Our result is slightly larger.   Fig. 25 and Fig. 26. The description of the plot is similar to Fig. 25. The left panel is for a = 0.06 fm and the right one for a = 0.04 fm.

Appendix C: Leading twist target mass correction
Unlike the light-cone ITD, the terms in the twist-2 OPE of the equal-time bilocal bilinear have trace terms, which are proportional to powers of hadron mass. At the level of twist-2 trace terms, such target mass effects have been calculated explicitly [65,66]. For the case of pion matrix element, such target mass corrected expressions are obtained from Eq. (24),Eq. (27) and Eq. (30) given by the replacement (P z z) n → (P z z) n n/2 k=0 (n − k)! k!(n − 2k)!
where n are even integer valued for the u − d pion PDF case. Including such correction terms in our analysis did not change the results (i.e., effectively the inferred values of valence pion moments from twist-2 OPE) well within their errors. However, there could be unaccounted target mass effects that originate from higher-twist terms and it is an expectation that their coefficients are O(1) or smaller, and hence suppressed as simple powers of m 2 π /P 2 z . We ignore any such effect in this computation.
Appendix D: Goodness of fits In Fig. 33, we plot the χ 2 /dof for our fits to the 2parameter and 4-parameter PDF ansatz, Eq. (46) in the main text. The description of the plot is similar to for F (x) = (n + α) log(x) + β log(1 − x) + log G(x). Now, we proceed towards doing a saddle point approximation in order to evaluate the leading term in the above integral in the limit of infinite n. The maximum of F (x) occurs at x = x 0 = 1 − β/n + O(1/n 2 ), which is less than 1 and hence within the domain of integration, and on the real axis. Thus, F (x) in the proximity of x = x 0 is (F3) Thus, the saddle point approximation gives the asymptotic dependence on n, from the first term in Eq. (F3) and an extra n from the change of variables in the last term of Eq. (F3) to perform the remaining Gaussian integral. The asymptotic series for x n in the limit of large-n is given by the standard multiplication correction factor which is a series in 1/n. FIG. 36. Plot shows the first three even moments obtained on a = 0.06 fm lattice in a model independent manner as described in Section VII. The description of the three plots is similar to that in Fig. 17, Fig. 18 and Fig. 19. In addition to the two values of P 0 z = 0, this plot includes P 0 z = 0 reference scale for the ratio. −(20) (39) with the same methodology as in Section VII.
In Fig. 37, we present results using fits to PDF ansatz. We refer the reader to Section VIII for explanations and methodology. In the top panel, the ratio matrix element with P 0 z = 0 is shown along with the bands resulting from fits to the 4-parameter ansatz. In the bottom panel, the best fit PDF using 4-parameter ansatz is shown (green) and compared to results using other non-zero P 0 z . The PDFs using different P 0 z remain more or less the same. The values of the exponent including results from P 0 z are α = −0.40  (46) . These results are to be compared with the entries in Table III. There is a tendency for P 0 z = 0 to pull the result for β higher; an opposite behavior wherein β increased when P 0 z is increased would have been more desirable. The inferred moments from model-dependent ansatz including P 0 z = 0 ratio with other ratio schemes give x v = 0.2470  (30) . The analysis of effective β that we discussed in Section IX gives β eff (n = 4) = 2.04 +(33)(50) (35) (50) that is again consistent with both β = 1 and β = 2.
Appendix I: Analysis imposing β = 2 in PDF ansatz In Section VIII, we used 2-and 4-parameter PDF ansatz in Eq. (46) to reconstruct the PDF that best describes the real-space lattice data. In that analysis, we kept β as a free parameter. Through that analysis, we found PDFs with β ≈ 1 or less to best describe the data. Here, we do the following; we take β = 2 as if it is a well established fact, and impose the constraint β = 2 in the 4-parameter ansatz and fit only α, s and t to minimize χ 2 . That is, even though there is a set of PDF that better describe the lattice data in the space of (α, β, s, t), we restrict now to a subspace (α, β = 2, s, t) and ask what PDFs within this subspace best describes the data.
Let us take a specific case of P 0 z = 0.43 GeV ratio ma- trix element on a = 0.06 fm lattice. We fit the lattice data in the range z ∈ [a, 0.5 fm]. The resultant χ 2 /dof for this 3-parameter ansatz was of course larger compared to the 4-parameter ansatz, but not very large either; for the 4-parameter ansatz for the same case, χ 2 /dof ≈ 25/36 while it was 56/37 for the 3-parameter one tried here. The resulting PDF is shown as the blue band in Fig. 38. For comparison, the unconstrained 4-parameter PDF result is also shown. The β = 2 result closely hugs the best-fit result and tries to lie within the 1-σ vicinity of it. The error bar on the constrained PDF is small because there only exists a small set of PDF with β = 2 that have a decent χ 2 . One might wonder if imposing β = 2 took our result closer to the ASV result [27]. For this, the ASV result is shown as the green dashed line Fig. 38. The β = 2 result actually misses the ASV result badly at intermediate x at the expense of agreeing well very close to x = 1. It is fascinating that this is actually due to robust tendency of the PDFs to pass through an approximate fixed-point at x = x * (≈ 0.6 specific to our data) that we discussed in Appendix G and determined by the first moment. Not surprisingly, the β = 2 fits resulted in values of first four moments as 0.254(5), 0.108(3), 0.057(2) and 0.034(2) for the case discussed, which compares well with the first four moments for the same case obtained using a full 4-parameter fit; namely, 0.245(8), 0.111(3), 0.064(4) and 0.040 (5). This analysis for fixed β = 2 assumes a very specific functional form for G(x) = 1 + s √ x + tx.
Thus, it is very much a possibility that by choosing some other flexible functional form for G(x), one might still be able to get β ≈ 2 and get better χ 2 . We do not explore this any further in this paper since effective β analysis addresses this in a better manner.