Pion Valence Structure from Ioffe-Time Parton Pseudodistribution Functions

We present a calculation of the pion valence quark distribution extracted using the formalism of reduced Ioffe time pseudo-distributions or more commonly known as pseudo-PDFs. Our calculation is carried out on two different 2+1 flavor QCD ensembles using the isotropic-clover fermion action, with lattice dimensions $24^3\times 64$ and $32^3\times 96$ at the lattice spacing of $a=0.127$ fm, and with the quark mass equivalent to a pion mass of $m_\pi \simeq 415$ MeV. We incorporate several combinations of smeared-point and smeared-smeared pion source-sink interpolation fields in obtaining the lattice QCD matrix elements using the summation method. After one-loop perturbative matching and combining the pseudo-distributions from these two ensembles, we extract the pion valence quark distribution using a phenomenological functional form motivated by the global fits of parton distribution functions. We also calculate the lowest four moments of the pion quark distribution through the"OPE without OPE". We present a qualitative comparison between our lattice QCD extraction of the pion valence quark distribution with that obtained from global fits and previous lattice QCD calculations.


I. INTRODUCTION
The key element of most predictions involving hard inclusive reactions in high-energy physics is the factorization theorem [1] of perturbative QCD. This factorization procedure separates the perturbatively calculable hardscattering quark and gluon dynamics from the nonperturbative bound-state dynamics, described by the parton distribution functions (PDFs) of the relevant hadrons.
For the nucleon, information about the quark PDFs can be obtained from the experimental data of deep inelastic scattering. Numerous experiments have been performed, and the functional form of the valence quark PDFs is well understood in various global fits [2][3][4][5][6].
On the other hand, the pion valence PDF has been extracted using the data from only a few pionic Drell-Yan experiments at CERN [7,8] and Fermilab [9]. The valence PDF of the pion is of particular theoretical interest, as the pion is the lightest QCD bound state and the Goldstone mode associated with the spontaneous breaking of chiral symmetry.
The accurate form of the pion's valence quark distribution therefore provides a testing ground for both QCD and QCD-based approaches in understanding the structure of hadrons. The experimental data of Refs. [7][8][9] have been analyzed in Refs. [10][11][12][13][14][15][16] to determine the pion valence distribution and compare, in particular, the results of these analyses at large x (fraction of the hadron's longitudinal momentum carried by the parton) with the predictions of QCD-based hard-gluon-exchange models [17][18][19].
Thus, the slope of the pion valence quark distribution as x → 1 may provide important information about the interplay of perturbative and nonperturbative aspects of quark dynamics in the valence region. Because of its importance, understanding the large-x behavior of the pion valence quark distribution is the goal of the approved experiment C12-15-006 at Jefferson Lab [27]. Intensive studies of the pion structure are also proposed for the future Electron-Ion Collider [28].
Another way of studying the behavior of the pion valence quark distribution may be provided by lattice QCD calculations which can be complementary to the global fit analyses of the cross-section data and can also serve as a discriminator between different model predictions.
To achieve such a goal, we need to go beyond the conventional moment calculations and obtain x-dependent parton distributions. Several new approaches that allow us to determine x-dependent parton distributions from lattice QCD have been proposed. These approaches include the path-integral formulation of the deep-inelastic scattering hadronic tensor [29], the inversion method [30], quasi-PDFs [31], good lattice cross sections [32,33], and reduced Ioffe-time pseudodistributions (or pseudo-PDFs) [34,35]. An analogous coordinate-space method has been earlier introduced for the calculation of light-cone distribution amplitudes [36].
Although significant achievements in the lattice QCD implementations of these approaches have been made in recent years [37][38][39][40][41][42][43][44][45][46][47][48], a proper understanding and controlling various sources of systematics in these calculations still require further exploration and theoretical development. The status of current lattice QCD calculations of the xdependent hadronic structure can be found in the following review articles [49][50][51]. Recently, an attempt to incorporate lattice QCD determination of PDFs together with experimental data to obtain the nonsinglet quark distribution of the nucleon has been discussed in [52].
Lattice calculations of the pion valence PDFs have been recently performed in Refs. [44,47,53] using the quasi-PDF [47,53] and the good lattice cross sections [44] approaches. In this paper, we present the first lattice calculation of the pion valence PDF using the approach based on reduced Ioffe-time pseudodistributions [34]. We discuss the limitations of our lattice QCD calculation of the pion valence quark distribution and compare our results with those in Refs. [44,47,53] and also those obtained from global fits as mentioned above.
The remainder of this article is organized as follows. In Sec. II, we briefly discuss the reduced Ioffe-time pseudodistributions approach and the necessity of the calculation in coordinate space. In Sec. III, we present numerical details of the calculation of hadronic matrix elements of the reduced Ioffe-time pseudodistribution to extract the pion valence quark distribution. We present the results of the extracted pion valence quark distribution in Sec. V and compare our result with different fits of the experimental data and other lattice calculations in Sec. VI. Finally, we summarize our results and outline the future directions of this method to obtain the pion valence quark distribution with controlled systematics.

II. BASICS OF THE IOFFE-TIME PSEUDO-DISTRIBUTIONS APPROACH
The unpolarized quark nonsinglet PDF is defined as a Fourier transform of the nonlocal matrix element taken on the light cone, e.g., for z ¼ ðz þ ¼ 0; z − ; 0 ⊥ Þ, with the momentum given by The combination ν ¼ p · z is called the Ioffe time [54], and Wðz; 0Þ is the gauge link in the fundamental representation. Its path goes along a straight line 0 → z. For general z, p and α, the Lorentz decomposition of this matrix element can be written as When z ¼ z − and α ¼ þ, the second function N ðν; z 2 Þ does not contribute; i.e., the twist-2 PDF is solely determined by the first function. On the lattice, we need to take a spacelike z. Choosing z ¼ z 3 and p ¼ ðE; 0 ⊥ ; p 3 ; Þ we can exclude the N ðν; z 2 Þ function and deal with the function Mðν; z 2 Þ that is called the Ioffe-time pseudodistribution (pseudo-ITD). The term "pseudo" [34] reflects the fact that one deals with the matrix element off the light cone, i.e., for nonzero z 2 . Another advantage of taking the time component of M α ðp; zÞ on the lattice is that such a choice also avoids the pitfall of renormalization constant mixing as described in Ref. [55]. Choosing a spacelike separation z brings in a serious complication of additional link-related ultraviolet (UV) divergences [56] that are absent when z is on the light cone. Fortunately, these divergences are multiplicatively renormalizable [57][58][59], i.e., form an overall factor Zðz 2 =a 2 Þ, where a is a UV regulator, such as the lattice spacing. In the quasi-PDF approach, such divergences are usually removed using various versions of the regularization-independent momentum subtraction (RI/MOM) method [60,62].
A different approach was proposed in [34]. One considers the reduced pseudo-ITD formed by the ratio of Mðν; z 2 Þ to its rest-frame p z ¼ 0 value Mð0; z 2 Þ. Since the UV factor Zðz 2 =a 2 Þ does not depend on ν, it disappears from the ratio Mðν; z 2 Þ. As a result, the latter is UV finite. Moreover, it is a renormalization group invariant quantity. Also, taking the ratio (3) removes not only the UV divergences, but also the part of the z 2 dependence associated with them.
Beside UV divergences, there are other sources of the z 2 dependence. In particular, the function Mðν; z 2 Þ contains higher-twist contributions related to the transverse-momentum distributions of quarks inside a hadron and reflected in the z 2 dependence of Mðν; z 2 Þ. For small z 2 , they appear as higher-twist contributions that are polynomial in Oðz 2 Λ 2 QCD Þ. Thus, one may expect that for small enough z 2 they may be neglected, and Mðν; z 2 Þ may be related to the light-cone Ioffe-time distribution [63] Qðν; μ 2 Þ. Such a relation is given by a factorization formula that involves just logarithmic lnðz 2 Þ dependence accompanied by a perturbatively calculable kernel.
An alternative proposal of taking such a ratio similar to that in Eq. (3) has been proposed in [61]. In that article, the authors claim that a vacuum matrix element of the same type of operator could be used in the denominator, instead of the rest frame hadron matrix element.
Furthermore, it is not unreasonable to suppose that the higher-twist z 2 dependence of Mðν; z 2 Þ and Mð0; z 2 Þ is similar, and the higher-twist impact on the ratio Mðν; z 2 Þ is much weaker than Mðν; z 2 Þ. In particular, if Mðν; z 2 Þ factorizes like Mðν; z 2 Þ ¼ MðνÞBðz 2 Þ, the ratio Mðν; z 2 Þ has no z 2 dependence. Unfortunately this idealized scenario is not true for QCD, but as was shown in [64], taking the ratio will always reduce the higher-twist contribution, in the limit ν → 0. An actual calculation [38], though performed in the quenched approximation, produced an almost z 2independent result for Mðν; z 2 Þ in the region z > 4a (i.e., in the region where one would expect higher-twist effects to be significant) and a logarithmic lnðz 2 Þ dependence in the region z ≤ 4a, where it was perfectly described by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [65][66][67] evolution.
Summarizing, the choice of Mðν; z 2 Þ as the basic object for lattice studies of parton distribution functions satisfies the following criteria: (i) Due to Lorentz invariance, the pseudo-ITD Mðν; z 2 Þ and, hence, the reduced pseudo-ITD Mðν; z 2 Þ depend only on the interval z 2 for a spacelike separation z between the quarks and the Ioffe time ν. Both z 2 and ν are Lorentz invariant. (ii) For all values of spacelike z 2 and for any contributing Feynman diagram, it can be shown [34] that the Fourier transform of Mðν; z 2 Þ and, hence, Mðν; z 2 Þ) with respect to ν has canonical support −1 ≤ x ≤ 1 in the variable x interpreted as the standard momentum fraction. (iii) The UV divergences associated with the gauge link are canceled in the reduced ITD, and the latter is a renormalization group invariant quantity. (iv) The short-distance behavior due to lnðz 2 3 M 2 Þ terms (where M is an infrared regulator and z ¼ z 3 ) present in Mðν; z 2 3 Þ and generating the perturbative evolution of parton densities is preserved in Mðν; z 2 3 Þ. In the z 2 3 → 0 limit, the reduced pseudo-ITD maps to the usual light-cone PDF and obeys the familiar perturbative DGLAP evolution with 1=z 2 serving as an evolution parameter. To get the matching condition between Mðν; z 2 3 Þ and the light-cone ITD Qðν; μ 2 Þ, one may use the operator product expansion (OPE) which is valid both for the numerator Mðν; z 2 3 Þ and the denominator Mð0; z 2 3 Þ. As discussed, the UV-singular factors present in these functions cancel together with the z 2 3 dependence associated with them. The remaining z 2 3 dependence corresponds to the DGLAP logarithms lnðz 2 3 M 2 Þ and higher-twist effects proportional to z 2 3 . In our particular case, we have an additional simplification that, in the local limit, the operator in Eq. (2) is a conserved vector current. As a result, the denominator does not bring an extra lnðz 2 3 M 2 Þ dependence. Eventually, Mðν; z 2 3 Þ is matched to the MS light-cone ITD by where Qðν; μÞ is the light-cone ITD whose Fourier transform with respect to ν gives the PDF fðx; μÞ at a factorization scale μ. The matching kernel Cðu; μ 2 z 2 Þ of the reduced pseudo-ITD to the MS ITD has been determined from one-loop calculations [59,[68][69][70] Cðu; where is the Altarelli-Parisi kernel [66] and The inverse of this formula will be needed for converting the lattice results for Mðν; z 2 Þ to the MS ITD Qðν; μÞ at a matching scale μ. Without loss of accuracy, this can be done by switching the reduced pseudo-ITD and the ITD and by changing the sign of α s .

III. NUMERICAL METHODS
The pion reduced pseudo-ITD is calculated on two lattice QCD ensembles with different physical volumes. These configurations were generated by the JLab/W&M Collaboration [71] using 2 þ 1 flavors of stout smeared clover Wilson fermions and a tree-level tadpole-improved Symanzik gauge action. The strange quark mass was set by requiring the ratio ð2M 2 K þ − M 2 π þ Þ=M Ω − to assume its physical value.
The fermion action includes one iteration of stout smearing with the weight for the staples given by ρ ¼ 0.125. This smearing procedure has the consequence that the tadpole-corrected tree-level clover coefficient c SW ¼ 1.2493 is very close to the nonperturbative value determined a posteriori with the Schrödinger functional method [71]. The ensemble parameters are listed in Table I. The lattice spacing of these ensembles, a ¼ 0.127 fm, was determined using the Wilson flow scale w 0 [72].
To ameliorate the contamination of excited states and improve the overlap of the interpolators onto boosted pions, we implement a combination of the Gaussian smearing [73] and momentum-smearing [74] techniques. The pseudo-ITD matrix elements are calculated using the summation method. By now, this method is well known in the lattice community. For completeness, some key points of the method are highlighted below. We refer the readers to [38,[75][76][77] for more details on the implementation of this method.
The summation method is related to the Feynman-Hellmann theorem. One considers a theory where the action is modified by the operator of interest By the Feynman-Hellmann theorem, a hadron matrix element of that operator can be found from a derivative of the energy of that hadron As was shown in [75], the derivative of the effective mass can be shown to be a ratio of correlation functions. This method has an advantage over other methods based around the Feynman-Hellmann method, which require generation of specialized configurations which scan λ. The derivative of the effective mass at λ ¼ 0 can be calculated using standard gauge configurations with λ ¼ 0.
The two-point and three-point correlation functions for a fixed pion momentum p and a fixed current insertion time t are written in terms of the standard pion interpolation field O Π : where T is the Euclidean time separation between the interpolating operators for the pion creation and annihilation operators, O Γ ðzÞ ¼ψð0ÞΓWð0; zÞψðzÞ, and we use Γ ¼ γ 4 . For a fixed z, summing over the current insertion time t, the matrix element is estimated from the large Euclidean time limit of the effective matrix element: where The leading excited-state effects can be parameterized by where Δ is the energy gap between the ground state and the lowest excited state. This method will also have significantly reduced excited-state contamination at large Euclidean separation compared to the typical ratio method, whose excited-state contamination decays as e −ΔT=2 . These decreased excited-state effects allow for particularly short time extents for matrix element extraction, as was demonstrated in [76]. For the summation method to be successful, many source-to-sink separations are required. The common sequential source technique would require a large number of propagator inversions for this to be practical. Instead, this calculation shall use a sequential operator to construct three-point correlation functions. The sequential operator H is defined as the solution to the system of equations, with suppressed spin and color indices X  I. The parameters for the JLab/W&M Collaboration ensembles used in this work: lattice spacing, pion mass, β, light and strange quark mass, spatial and temporal size, and the numbers of configurations. The a127m415L ensemble contains ten independent streams of 256 configurations each while the a127m415 ensemble contains a single stream.
calculated by replacing a single propagator in a standard two-point function calculation with the sequential operator. As mentioned earlier, this calculation utilizes the momentum-smearing procedure [74] to improve the signal of the moving states. Three values of the momentumsmearing parameter ζ are used, including 0. For each of those momentum-smearing parameters, two motifs of smearing are used. The source interpolating field is always smeared, but the sink quarks will be either smeared, called smeared-smeared (SS), or left as points, called smearedpoint (SP). Typically the SS correlation functions will have lower excited-state contamination. On the other hand, the SP correlation functions will typically have less statistical noise.
A fit of the data to Eq. (14) is used for each correlation function, holding the ground-state matrix element and effective energy gap fixed between each correlation function of the same p and z. Specifically, a fit is performed on N different effective bare matrix elements with different smearing setups, M eff j , to the form with (2N þ 2) fit parameters where j labels the N different smearings. The fit parameters will be chosen with a weighted χ 2 minimization which employs the full covariance matrix. The different momentum-smearing parameters only improve the signal-to-noise ratio for a certain range of momenta states and decrease the signal-to-noise ratio for other momenta. The T range as well as which of these correlation functions are used in the fit are varied to minimize the χ 2 per degree of freedom (χ 2 =d:o:f:). All statistical errors for the correlation functions and matrix elements are estimated using the jackknife resampling technique. Examples of these fits for the pion matrix elements are plotted in Figs. 1(a) and 1(b). The bare matrix elements for the rest frame are shown in Fig. 2. A unique feature of pion correlation functions at zero momentum is a constant signal-to-noise ratio in T allowing these matrix elements to be significantly more precise than they would be for other hadrons. The value of these matrix elements decays exponentially in z=a. This feature is generated by the renormalization of the Wilson line. In perturbation theory, this behavior appears as a power divergence for small z=a. This exponential behavior will appear in matrix elements for all momenta, but without the constant signal-to-noise ratio. The matrix elements for large distance become very small and with the exponentially growing statistical error at p ≠ 0, they can be increasingly difficult to resolve.
In our calculation, the largest momentum along the z direction we use for both of the ensembles is p max ¼ 3ð2π=LaÞ. We use z max ¼ 6a (0.76 fm) and 8a (1 fm) for the a127m415 and a127m415L ensembles, respectively. Even with these relatively large separations, there does not appear a noticeable sign of higher-twist effects within the limitation of present statistics. Consequently, when calculating the moments of the PDF or of the MS ITD, both of which can be calculated for each z independently, there does not appear to be any significant dependence on the value of z used as will be demonstrated in Sec. IV. If one attempts to determine the PDF from the ITD using a limited set of data, i.e., z ≤ 4a, or using the full set of data, then the results will be consistent with each other but with larger variance for the smaller set of data. This feature is particularly apparent in the low-x region due to the shortened Ioffe-time extent of the ITD. It is important to note that the size of potential highertwist effects must be confirmed before trusting results at any separation. As is frequently done, labeling higher-twist effects as Oðz 2 Λ 2 QCD Þ for pseudo-ITDs or OðΛ 2 QCD =p 2 z Þ for quasi-PDFs only estimates the size of these effects and the data must be checked for the presence or lack of these effects before results can be trusted. For the largest separation in this analysis, we have z 2 Λ 2 QCD ∼ 1, but the reduced pseudo-ITD appears to have successfully removed the higher-twist effects in this range of Ioffe time.
The reduced-ITDs for both ensembles before any perturbative matching are shown in Fig. 3. One can see that the most data group around some ν-dependent curve, with a rather small scatter. Excluded from this pattern are the p z ¼ 3 data points on the a127m415 ensemble, which are visibly outside from the other data points that follow a somewhat regular distribution in ν in Fig. 3. It is worth noting that while p z ¼ 3 corresponds to a momentum of ∼0.91 GeV in physical units on the a127m415L ensemble, it corresponds to a momentum of ∼1.22 GeV on the a127m415 ensemble with smaller volume. In fact, these data points are still within less than ∼2σ away of the other data points, and their inclusion does not affect the subsequent result of the fit to extract the pion valence PDF.

IV. MOMENTS OF THE PION PDF
As was described in [64], the reduced pseudo-ITD can be used to calculate the moments of the PDF. By comparing the Taylor expansions with respect to ν of Eq. (4), one can derive a multiplicative relationship between the moments of the pseudo-PDF, b n ðz 2 Þ, and the moments of the MS PDF, a n ðμ 2 Þ: where C n are the Mellin moments of the matching kernel Cðu; μ 2 z 2 Þ with respect to u. To next-to-leading-order (NLO) accuracy, the moments are given by where are the well-known moments of the Altarelli-Parisi kernel and By completely avoiding the inverse problem [64], this procedure allows for an understanding of the PDF's structure before any potential systematic errors arising from the matching convolution and Fourier transforms are incurred. In principle, this method, coined as "OPE without OPE" [78], can be used to determine any moment of the PDF in sharp contrast to the traditional method which is based on local matrix elements. The latter is limited by FIG. 2. The bare matrix element calculated for the case of p ¼ 0. The constant signal-to-noise ratio in T allows for these pion matrix elements to be extremely precise compared to other hadrons. The exponential decay in z is a feature caused by the renormalization of the Wilson line.
FIG. 3. Real component of the reduced pseudo-ITD obtained from ensembles a127m415 and a127m415L for z max ¼ 6a and 8a, respectively. The largest momentum used for both the ensembles is p max ¼ 3ð2π=LaÞ. The triangle-up (△) symbols indicate the reduced pseudo-ITD matrix elements Mðν; z 2 Þ extracted from the a127m415 ensemble and the hexagon (⬡) symbols denote those for the a127m415L ensemble. These data points represent the reduced pseudo-ITD before any perturbative matching is performed. the appearance of power divergent mixing due to the reduced rotational symmetry of the lattice as well as issues related to the signal-to-noise ratio. In practice, however, only the lower moments will have a resolvable signal. Larger Ioffe-time extents and a finer resolution in Ioffe time, both of which require finer lattice spacings, are necessary before higher moments can be obtained.
The lowest four moments will be extracted from the a127m415L ensemble data by inverting the Vandermonde matrix as was performed in [64,77]. The imaginary and real components are used individually to calculate the odd and even moments, respectively. The imaginary component of the reduced pseudo-ITD obtained from the a127m415L ensemble is shown in Fig. 4. Due to the larger uncertainty of the p z ¼ 3 data points from the a127m415 ensemble and the deviation from the Ioffe-time distribution of the other data points as shown in Fig. 3, this ensemble does not allow us to extract moments in a reliable way. Therefore, we extract moments using the real and imaginary components of the reduced-ITD only from the a127m415L ensemble in what follows.
The results for these moments, determined from data with different z independently, are shown in Figs. 5-8, each one called "pseudo-PDF moment" before matching and "PDF moment" after matching. For the matching relationships, we choose μ ¼ 2 GeV and the value of α s ð2 GeVÞ ¼ 0.303. This value of the coupling is taken from the evolution used by the LHAPDF [79] for the dataset cj15nlo from the CTEQ-Jefferson Lab Collaboration [80]. The imaginary component of the pseudo-ITD calculated with the lowest two z values has a completely linear behavior within their short Ioffe-time range. This leads to an inaccurate determination of the third moment, which also slightly affected the calculation of the first moment. The first moment for these two separations is instead calculated with a linear fit. Similarly, the real component of the reduced pseudo-ITD calculated with the lowest three z values only exhibits a quadratic behavior and has an illconstrained value for the fourth moment. In this case, those results did not appear to affect the quality of the second moment determination, so the results are kept for the second moment. For the third and fourth moments, these poorly constrained results are dropped from the following analysis. FIG. 4. Imaginary component of the reduced pseudo-ITD obtained from the ensemble a127m415L for z max ¼ 8a. The largest momentum is p max ¼ 3ð2π=LaÞ.
FIG. 5. The first moments of the pion pseudo-PDF and of the PDF calculated from the a127m415L ensemble are shown. The moments are shown as functions of the original separation ðz=aÞ 2 of the reduced pseudo-ITD data used to calculate them with small offsets for better visibility. After the matching procedure, the dependence on the separation is significantly reduced. The largest deviation occurs for the lowest separation, which is most susceptible to discretization errors.
FIG. 6. The second moments of the pion pseudo-PDF and of the PDF calculated from the a127m415L ensemble are shown. The moments are shown as functions of the original separation ðz=aÞ 2 of the reduced pseudo-ITD data used to calculate them with small offsets for better visibility. After the matching procedure, the dependence on the separation vanishes within statistical precision, showing the lack of higher-twist effects in this moment at this level of precision. The green band represents the PDF moment from a weighted average. After the matching procedure in Eq. (16) has been applied, there is no significant large-z dependence on the separation z from which the data originated for the second, third, and fourth moments. On the other hand, a z 2 dependence does appear for the first moment calculation, particularly between the z ¼ 1a and 1a < z ≤ 8a calculations. This moment at z ¼ 1a point is the most sensitive to the lattice spacing errors. In addition, this point was determined with a simple linear fit and therefore it can be affected by different systematics.
For the first moment the residual z dependence at larger ðz=aÞ 2 , which may be due to higher-twist effects, can be modeled. These effects are either polynomial (i.e., z 2n ) or of the form z 2n lnðz 2 μ 2 e 2γ E þ1 4 Þ. The logarithmic terms arise from logarithms in the matching coefficients C n ðz 2 μ 2 Þ when applied to pseudo-PDF moments that contain higher-twist effects. As we can see these terms are suppressed by α s and they are expected to be smaller than the simple polynomial terms. Due to its large potential discretization errors, the z ¼ a data point is neglected in this fit. The moments are fit to four different functional forms: ; and the results of these fits are shown in Table II and are plotted in Fig. 9. The value of Λ QCD ¼ 300 MeV was used, but this choice is made solely for the magnitude of the coefficients of the Oðz 2 Λ 2 QCD Þ terms to be estimated. From the value of the coefficients in Table II, the higher-twist effects in this moment are an order of magnitude smaller than what a naïve Oðz 2 Λ 2 QCD Þ estimate would have suggested. The first moment of the PDF dominates the low-ν behavior of the pseudo-ITD and ITD, which is precisely the region where the reduced pseudo-ITD reduces higher-twist effects. The systematic error introduced by these highertwist effects would be comparable to, or smaller than, the other systematic errors of finite lattice spacing and unphysical pion mass. If there had appeared a sufficiently strong polynomial z 2 behavior in enough moments, one could use the inverse Mellin transform of that behavior to estimate the higher-twist effects in the reduced pseudo-ITD. With only four moments and only one of them showing any discernible higher-twist effects, this inverse transform will not be reliable. This lack of large higher-twist effects in the The moments are shown as functions of the original separation ðz=aÞ 2 of the reduced pseudo-ITD data used to calculate them with small offsets for better visibility. The data from the lowest three z values in this calculation only show linear effects due to the short range of Ioffe time they span. As a result, they do not have any signal for the fourth moment and are not shown here. After the matching procedure, the dependence on the separation vanishes within statistical precision, showing the lack of large higher-twist effects at this level of precision. The green band represents the PDF moment from a weighted average. moments, even with separations as large as z ¼ 8a, justifies the use of these data for extracting the PDF.
We summarize our calculation of the moments of the pion PDF in Table III. For the moments which do not show signs of higher-twist effects, we take a covariance weighted average of the results for each separation. For the first moment, we state the value from the fit a latt 4 with the highertwist contamination removed. The fit a latt 4l has large values of the variance of the fit parameters, a small number of degrees of freedom, and the oscillatory nature of the final result. We believe this result obtained from a latt 4l is overfit. We refer the readers to the previous calculations [81][82][83][84][85][86][87] of such moments for a comparison with the calculated moments in this work.
We shall calculate the moments from the extracted pion valence PDF, in the next section, and compare with those obtained from in the NLO QCD analysis [14] from the Fermilab E-615 pionic Drell-Yan data [9] at a scale of 5.2 GeV. Note that the odd moments listed in Table III, which are extracted from the imaginary component of the reduced pseudo-ITD, are related to q v ðxÞ þ 2qðxÞ distribution as shown in Ref. [38]. Therefore, our results for the odd moments are not directly comparable with that calculation.

V. EXTRACTION OF THE PION VALENCE DISTRIBUTION
It can be seen from Fig. 3 that the z 2 dependence indeed cancels out to a great extent in the reduced pseudo-ITD of Eq. (3) without spoiling the ν dependence in the Ioffe-time distribution which governs the shape of the pion PDF. As discussed earlier, for small z 2 the function Mðν; z 2 Þ should contain lnðz 2 Þ singularities related to the perturbative evolution of the PDFs. Such a logarithmic z 2 dependence is clearly seen in the z ≤ 4a data of Refs. [38,39], performed in the quenched approximation, though at a finer lattice spacing a ≃ 0.093 fm. For large z ≥ 6a values, they found the data in practice did not depend on z. Thus, one can explicitly identify the region z ≤ 4a, where one may rely on the perturbative evolution. As demonstrated in [69], the matching procedure applied in this region to the points for Mðν; z 2 Þ converts their lnðz 2 Þ dependence into the lnðμ 2 Þ dependence of the MS light-cone ITDs Qðν; μÞ on the MS scheme subtraction parameter μ 2 . Just like in Ref. [77], for the matching of the reduced pseudo-ITD to the MS ITD at a particular scale μ, we perform an inversion of Eq. (5) simply by switching Qðν; μÞ and Mðν; z 2 Þ and changing the sign of α s . This gives In other words, applying the matching procedure for an appropriate value of α s , the z 2 dependence of the original small-z 2 data for Mðν; z 2 Þ should be compensated by the FIG. 9. The first moment of the PDF calculated from the a127m415L ensemble is shown. The different bands correspond to different models of the higher-twist effects. These effects are significantly smaller than what an Oðz 2 Λ 2 QCD Þ estimation would have provided. The value hx 1 i ¼ 0.2541ð26Þ denoted by the green band in a latt 4 fit is quoted in Table III. ln z 2 term, and one should get practically z 2 -independent data points for Qðν; μÞ.
In the present calculation, due to significantly larger uncertainties compared to those in [38], no systematic logarithmic dependence of the data on z 2 is visible. Therefore, it is not possible to determine what is the length scale z 0 below which one may rely on perturbative evolution of the reduced pseudo-ITD and what is the value of α s associated with this evolution. As guidance for a particular choice of α s one may use the idea that using the correct choice of α s in the matching formula (20) should produce the least scatter for the small-z 2 points from a universal curve in ν, at least up to discretization effects.
The convolution is performed on a polynomial fit of the data for Mðuν; z 2 Þ, and Eq. (20) is applied for each z independently. We choose μ ¼ 2 GeV, and the value of α s ð2 GeVÞ ¼ 0.303 has been taken from the evolution used in Ref. [79] as mentioned earlier. It important to note that, in the reduced-ITD approach, the relevant scale for converting the matrix element to the MS scheme is the separation between the quark fields z 2 , not the hadron's momentum p z . To investigate the systematics of the oneloop matching on our choice of α s at a particular scale, we vary α s by 10% and estimate its effect as a source of systematic uncertainty, as shown in Fig. 10. One can show that only the data points at large ν have as large as 5% change in their central values and are still statistically consistent between different choices of α s . Other points at low Ioffe time have less than 1% change for the variation in α s ¼ 0.303 AE 0.030. The small differences in the matched reduced-ITDs of the two ensembles originating from 10% change in α s are propagated as systematic uncertainty in the subsequent analyses. The matched ITDs for the two ensembles at the matching scale of 2 GeV are shown in Fig. 10. The matched data points in Fig. 10 are also seen to be located along a single curve that is a function of ν, with rather small fluctuations. Again, the exceptions are the p z ¼ 3 data points from the a127m415 ensemble.
Having data from two ensembles that have different volumes, one may wish to use them to study effects due to finite volume. We do not see any finite volume effect in the reduced-ITDs with the exception of p z ¼ 3 data obtained from the a127m415 ensemble which is not trustworthy for the small volume. Moreover, the present case of only two volumes is not an ideal situation for an infinite-volume extrapolation. In our particular case, there are few pairs of data from the two ensembles that correspond to the same Ioffe time ν.
To extract a single PDF combining the results from the two ensembles, we will perform a simultaneous and correlated fit to these two datasets. The number of configurations of the two ensembles we use are different. Therefore, to perform a simultaneous fit of the matched Ioffe-time distributions, equal numbers of bootstrap samples are generated from the two ensembles. Because of the reason mentioned earlier, we exclude the p z ¼ 3 data from the a127m415 ensemble in this fit. Because the functional form of the ITD QðνÞ is not known a priori, we implement a "z expansion" [88,89] fit to the datasets, reflecting the analyticity of QðνÞ. We also investigate whether there is residual z 2 dependence in the matched Qðν; μÞ distribution. This is achieved by adding z 2 -dependent terms in the fit function. Finally, a term which describes a potential volume dependence is added. There are no model calculations for the form of the finite volume corrections to the matrix element in Eq. (1), unlike the two current matrix element [90]. Instead, we take a simple exponential volume dependence where the relevant distance is the difference between the lattice size L and the length of the Wilson line z. Namely, we use the following form: where Note that unlike the form factors we do not have cuts in the complex plane for the calculated reduced pseudo-ITD and we simply choose a dimensionless number ν cut ¼ 1.0 in Eq. (21). In fact, other choices are possible with the final results being unaffected.
One can readily see from Eq. (22) that, for fixed values of λ k , a larger ν cut dictates a slower falloff of the distribution. However, if one allows to vary λ k in the fit along with different choices of ν cut , the fit parameters are also changed accordingly such that the red band shown in Fig. 10 remains unchanged. The value λ 0 ¼ 1.0 in the fit is fixed by the normalization in Eq. (3) at ν ¼ 0 and does not change by the perturbative matching. We limit the k max in our fit to the value in which additional term in the z expansion has no effects on the fit and obtain k max ¼ 4. The fit parameters are listed in Table IV. The smallness of the fit parameters c 1 and c 2 reflect the fact that the residual z 2 dependence is negligible as discussed earlier. The red band in Fig. 10 represents the ITD in the limit of infinite volume and vanishing z 2 contribution. The outer red band indicates the systematic uncertainty in the z expansion fit of the reduced-ITD introduced by 10% variation in α s at μ ¼ 2.0 GeV.
We now use the Qðν; μ ¼ 2 GeVÞ ITD from the above fit to extract the pion valence quark distribution. By definition, QðνÞ and the valence quark distribution of the pion q π v ðxÞ are related by and the quark distribution is given by the inverse Fourier transform Therefore to convert QðνÞ into a function of x, one should, in principle, know QðνÞ for all ν. In our lattice QCD calculation, we are restricted to ν max ∼ 6 for a discrete set of integer ν. Thus, the extraction of the PDF using Eq. (24) from lattice-calculated data constitutes an ill-posed inverse problem. To our knowledge, a reliable direct inverse Fourier transform to extract the PDF using lattice QCD data is currently a formidable task. An important constraint serving as additional information is that the valence distributions of the nucleon and the pion are smooth functions of the momentum fraction x in the region 0 < x < 1, with support only in that region, and the pseudo-ITD is related by a cosine transform to the pion valence quark distribution [38,63] QðνÞ ¼ In the spirit of the functional forms used in global fits of PDFs, we insert into Eq. (25) and numerically perform the integration, where N is the normalization such that We use the numerical fitting program ROOT [91] to fit bootstrap samples of the matched pseudo-ITD to the q π v ðxÞ distribution. We find that the ρ ffiffi ffi x p or the γx term has no effect in the fit and we do not obtain any signal of the fit parameters ρ and γ. Therefore we drop these terms and adopt in our calculations the following simple functional form for the PDF: where the beta functions in the denominator ensure that the normalization condition in Eq. (27) is met. We obtain the fit parameters α ¼ −0.48ð14Þ stat ð4Þ sys ; with the χ 2 =d:o:f. about 1.9. For example, the inclusion of the γx term in the PDF fit yields the following fit parameters and therefore the extracted PDF remains essentially unchanged: In Eqs. (29) and (30), the numbers in the first uncertainty in the parentheses of the fit parameters are statistical uncertainties and the second is obtained from fitting the matched Ioffe-time distribution with the 10% variation in the value of α s . We present the extracted PDF q π v ðxÞ from this fit in Fig. 11(a) and the xq π v ðxÞ distribution in Fig. 11(b). The fit to the data returns a well-constrained value, α ¼ −0.48ð14Þ, for the small-x behavior of the PDF corresponding to the slope of the relevant Regge trajectory. However, we stress that, with the present resources of lattice QCD calculations, a precise and accurate determination of the low-x behavior of the PDFs is not accessible. Specifically, one can argue that the Fourier transform at the Ioffe time ν is related to the region around the inverse of the Bjorken variable x B , i.e., ν ¼ 1=x B [33]. Therefore, to obtain a reliable estimate of the low-x behavior of the PDFs, one requires knowledge of the ITD at large ν. The PDF model parameters of the fit are highly correlated, as is evident from Figs. 11(a) and 11(b). The uncertainty at x ¼ 0.65 shrinks significantly due to the correlation between the fit parameters. This feature of shrinking uncertainty at different x values has also been observed in the calculation of nucleon PDFs using pseudo-PDF approach [77]. It is a feature of these highly correlated fits to have regions with small statistical errors.
For a comparison with the OPE without OPE calculation of the moments in Table III, we can take the Mellin transformation of this PDF result described by the fit parameters in Eq. (29). At the scale of μ ¼ 2 GeV, the first four moments are hxi ¼ 0.188ð56Þ, hx 2 i ¼ 0.081ð29Þ, hx 3 i ¼ 0.046ð19Þ, and hx 4 i ¼ 0.030ð14Þ. Note that any discrepancy between these numbers and those in Table III may arise from a number of different reasons. Firstly, the odd moments quoted in Table III contain small contribution from the antiquarks. Secondly, the dataset used in the PDF extraction includes data from both ensembles while those in Table III only include data from the a127m415L ensemble. We also present a comparison between the moments extracted from our pion valence PDF fit with those extracted in the NLO QCD analysis [14] using the Fermilab experimental data [9] at a scale of 5.2 GeV in Table V. We note that this lattice QCD calculation is performed at an unphysical pion mass of m π ∼ 415 MeV. A proper comparison with the QCD analysis of the experimental data and lattice QCD calculation can be made when a continuum and infinite volume extrapolation to the lattice data near the physical pion mass is performed.

VI. COMPARISON WITH OTHER DETERMINATIONS
This lattice QCD calculation using the pseudo-ITD approach is performed at a relatively heavy pion mass (m π ≃ 415 MeV) and on a relatively coarse lattice spacing of a ¼ 0.127 fm. Repeating similar calculations on several other lattice ensembles to determine the pion mass TABLE V. Comparison between the moments of the PDF extracted from the pion valence distribution in this work and those obtained in the NLO QCD analysis [14] using the Fermilab experimental data [9] at a scale of 5.2 GeV. shows the xq π v ðxÞ distribution. The initial scale for evolving the PDF to a higher scale is μ ¼ 2 GeV, as described in Sec. VI. For the fits of PDFs, we use the covariance matrix obtained from the z-expansion fit to generate bootstrap samples in the Ioffe-time range of 0 < ν < 4.71 for which lattice QCD data points exist. We do not perform any extrapolation to the data points outside this region of the Ioffe time using the results of the z expansion. We have checked that there is no dependence on the number of bootstrap samples in our fit. The statistical uncertainty band is obtained from the fits to the bootstrap samples of the data. The inner red band is obtained from a simultaneous fit to the matched ITDs at μ ¼ 2.0 GeV on these two ensembles in the limit of infinite volume. The outer red uncertainty band in the extraction of PDFs is obtained as a source of systematic uncertainty by calculating the difference between the simultaneous fit to the matched ITDs for α s ¼ 0.303 AE 0.030. dependence, quantify the severity of discretization errors and provide a more precise estimation of the effect of finite volume in order to obtain the pion valence PDF in the continuum limit is under way. However, with the present calculation, we proceed with a qualitative comparison of various global fits of the pion PDFs and three previous lattice QCD determinations.
For a comparison with q π v ðxÞ determined from the Drell-Yan experimental data in Ref. [9], we evolve our determination of the pion PDF to an evolution scale of μ 2 ¼ 27 GeV 2 starting from an initial scale of μ 2 0 ¼ 4 GeV 2 as shown in Figs. 11(a) and 11(b). As expected, the evolution to a higher scale shifts the peak of the xq π v ðxÞ distribution toward smaller values of x and a more convex-up behavior of the distribution as x → 1 is seen compared to the xq π v ðxÞ at the initial scale of our calculation.
We see a quantitative agreement of our extracted pion PDF with the NLO global fits in [14,16] to the Drell-Yan experimental data in [9] in the x ≳ 0.7 momentum fraction region as presented in Fig. 12. It should be noted that our determination of q π v ðxÞ has a distinctive deviation from the predictions of QCD-based hard-gluon-exchange perturbative models [17][18][19], which predicts a faster ð1 − xÞ 2 falloff of the pion PDF at large x. Such a ð1 − xÞ 2 falloff at large x was also obtained in the experimental data analysis in Ref. [15], where the authors included next-to-leadinglogarithmic threshold soft-gluon resummation effects in the calculation of the Drell-Yan cross section. A good way to visualize the discrepancy in the large-x region between the pion PDF extracted in our calculation with the experimental data and various global fits can be demonstrated by plotting xq π v ðxÞ as a function of x. We present such a plot in Fig. 13.
We now compare the determination of the pion valence quark PDF presented in this paper with the previous three lattice QCD determinations using the quasi-PDF approach [47,53] and the lattice cross-sections approach [44].
In [47], a careful and systematic investigation was performed. In that paper, the matrix element was renormalized using the RI/MOM-inspired approach [60,62], in which the hadron matrix element is divided by a quark matrix element of the same operator. Thus, just like in the reduced pseudo-ITD method, one deals with a ratio of the original matrix element and another matrix element that has the same ultraviolet divergences.
However, the nonperturbative z 2 behavior of the denominator factor in these two approaches is different. In particular, according to the analysis in Ref. [47], the RI/ MOM factor Z RI=MOM ðzÞ in the z ≳ 0.5 fm region shows a formation of a constituent quark mass m scr ∼ 300 MeV leading to a suppression of Z RI=MOM ðzÞ by an extra e −m scr jzj factor.
The pion rest-frame matrix element Mð0; z 2 Þ used in our calculation has a much faster decrease with jzj than Z RI=MOM ðzÞ of Ref. [47]. Numerically, it is very close to the nucleon rest-frame matrix element of Ref. [77]. As argued in Ref. [34], the fast falloff of Mð0; z 2 Þ reflects the finite size of the relevant hadron, i.e., the nonperturbative effects related to quark confinement. In the OPE language, FIG. 12. Comparison of the pion q π v ðxÞ distribution with the leading-order (LO) extraction from Drell-Yan data [9] (gray data points with uncertainties), next-to-leading-order (NLO) fits [14][15][16] (green band, maroon curve, and blue band). This lattice QCD calculation of q π v ðxÞ is evolved from an initial scale μ 2 ¼ 4 GeV 2 at NLO. All the results are evolved to an evolution scale of μ 2 ¼ 27 GeV 2 . The outer red uncertainty band shown in the q π v ðxÞ distribution is obtained from the variation in the choice of α s during the one-loop perturbative matching as described in Sec. V.
FIG. 13. Comparison of the pion xq π v ðxÞ distribution with the LO extraction from Drell-Yan data [9] (gray data points with uncertainties), NLO fits [14][15][16] (green band, maroon curve, and blue band). This lattice QCD calculation of q π v ðxÞ is evolved from an initial scale μ 2 0 ¼ 4 GeV 2 at NLO. All the results are evolved to an evolution scale of μ 2 ¼ 27 GeV 2 . The outer red band shown in the q π v ðxÞ distribution is obtained from the variation in the choice of α s during the one-loop perturbative matching as described in Sec. V and by calculating the variation in the fitting of the matched Ioffe-time data using the PDF parametrization in Eq. (28). the associated z 2 dependence corresponds to higher-twist contributions. As pointed out in Sec. II, one of the aims of using the reduced ITD is to cancel unwanted higher-twist effects.
We note that the use of the reduced ITD corresponds to a gauge-invariant renormalization prescription that avoids the pathological systematics of fixed gauge renormalization. In fact, previous calculations of quasi-PDFs [53,92] have shown slight discrepancies which depend on the renormalization scheme used for the matrix element and intermediate schemes used in the matching relationships. The systematic errors introduced by these choices can be avoided by calculating a renormalization group invariant quantity Mðν; z 2 Þ and then matching it to the MS PDF.
Another difference among the different lattice calculations is the treatment of the inverse problem. Our point is that inverse problems can only be solved by adding additional information. The quasi-PDF calculation in [47] adds this information in a way analogous to the present calculation and the lattice cross-sections calculation in [44]. The main idea is to parameterize the PDF in terms of a few model-dependent parameters and then fit the position-space matrix element using that functional form.
The quasi-PDF calculation performed in Ref. [53] instead attempts to directly perform the inverse Fourier transform. To remove unphysical oscillations caused by the ill-posed inverse, they use a "derivative method" [93]. As was shown in [64,94], the derivative method does not alleviate the ill-posed inverse problem especially with such short Ioffe-time ranges, as only a few nonzero points exist in that quasi-PDF calculation [53].
Of notable interest, our calculation of xq π v ðxÞ, shown in Figs. 11 and 13, illustrates a peak of the distribution in a region x < 0.40. This is consistent with all the global analyses of the pion valence distribution, wherein xq π v ðxÞ is peaked below x ¼ 0.40. This feature also occurs in the lattice cross-sections calculation in [44] and the quasi-PDF calculation in [47]. On the other hand, the quasi-PDF calculation using the derivative method performed in [53] peaks at a somewhat larger value of around x ¼ 0. 50.
In comparison with the determination of q π v ðxÞ in Ref. [44], which was performed on the same ensemble as the larger of our lattices, we see the large-x behavior of the distributions is in agreement within their uncertainties. In particular, q π v ðxÞ extracted in [44] gives β ¼ 1.93ð68Þ, and the value β ¼ 1.08ð41Þ is obtained in the present calculation. We note, however, that the choice of the initial scale of 1 GeV performed in Ref. [44] was rather arbitrary because of the absence of the NLO perturbative matching. The difference in q π v ðxÞ obtained in these two calculations based on short-distance factorization remains to be investigated.
In Fig. 14, we present a comparison between the lattice QCD extractions of pion valence quark distribution made in Refs. [44,47,53] with our calculation that uses the pseudo-ITD approach. For the calculation in Ref. [44], the PDF is evolved to μ ¼ 4 GeV assuming an initial scale of 2 GeV in Fig. 14. The PDF in Ref. [47] is calculated using pion momentum p z ¼ 1.29 GeV, the RI/MOM scale is fixed at 1.93 GeV and the PDF is estimated at μ ¼ 3.2 GeV. To illustrate the difference between all the calculations in the x ≳ 0.4 region, we also present the xq π v ðxÞ distributions from these lattice calculations in Fig. 14(b).

VII. SUMMARY AND OUTLOOK
In this paper, we have presented the first lattice QCD calculation of the pion valence distribution using the Ioffetime pseudodistributions approach. In our calculation, we have used combined results from two different ensembles, with the same lattice spacing a ¼ 0.127 fm but with different lattice sizes of 24 3 × 64 and 32 3 × 96. We have also performed the first lattice QCD calculation of the fourth moment of the pion PDF, which was previously inaccessible from local matrix elements.
FIG. 14. Comparison of the pion valence q π v ðxÞ and xq π v ðxÞ distributions estimated near μ ¼ 4 GeV scale using quasi-PDFs in Refs. [47,53] and good lattice cross sections in Ref. [44] with the pseudo-ITD approach in this calculation.
We have performed the one-loop perturbative matching of the pseudo-ITDs to light-cone ITDs at the scale μ ¼ 2 GeV. Then we have combined the pseudodistributions from the two ensembles using a model-independent scheme.
To extract the pion valence quark distribution, we have assumed a functional form motivated by those used in phenomenological global fits of parton distribution functions and obtained the parameters of this form using our data for the matched light-cone ITD. This approach allows us to avoid solving the ill-posed problem of the inverse Fourier transform from the ITD to PDF. We made a qualitative comparison between our lattice QCD extraction of the pion valence quark distribution with those obtained from global fits and previous lattice QCD calculations.
It should be noted that, in our lattice calculation, the region z ≲ 0.5 fm where one can rely on the perturbative matching corresponds just to four lattice separations at this lattice spacing. Therefore, a calculation with finer lattice spacing is warranted and is the goal of future efforts. On the other hand, our analysis (even keeping in mind the limitations of the present statistics) shows no evidence for higher-twist effects in the reduced pseudo-ITD up to z ¼ 1 fm, which we took as a justification to use the data for extracting the PDF up to z ¼ 8a on the larger ensemble.
A natural extension of our investigation would be a calculation on several other ensembles with different pion masses, lattice spacings, and volumes. We expect that the results of such upgraded lattice QCD calculations will eventually become important additional input in global analyses to obtain a precise knowledge of the large-x behavior of the pion valence quark distribution.