First direct lattice-QCD calculation of the $x$-dependence of the pion parton distribution function

We present the first direct lattice-QCD calculation of the Bjorken-$x$ dependence of the valence quark distribution of the pion. Using large-momentum effective theory (LaMET), we calculate the boosted pion state with long Wilson link operators. After implementing the one-loop matching and meson mass corrections, our result at $m_\pi \approx 310$ MeV is in agreement with those extracted from experimental data as well as from Dyson-Schwinger equation in small $x$ region, but a sizeable discrepancy in the large $x$ region. This discrepancy provides a nice opportunity to systematically study and disentangle the artifacts in the LaMET approach, which will eventually help to discern various existing analyses in the literature.

We present the first direct lattice-QCD calculation of the Bjorken-x dependence of the valence quark distribution of the pion. Using large-momentum effective theory (LaMET), we calculate the boosted pion state with long Wilson link operators. After implementing the one-loop matching and meson mass corrections, our result at mπ ≈ 310 MeV is in agreement with those extracted from experimental data as well as from Dyson-Schwinger equation in small x region, but a sizeable discrepancy in the large x region. This discrepancy provides a nice opportunity to systematically study and disentangle the artifacts in the LaMET approach, which will eventually help to discern various existing analyses in the literature.

Introduction:
The pion plays a fundamental role in QCD. As the lightest meson and the Goldstone boson associated with dynamical chiral symmetry breaking, it provides an important testing ground for our understanding of nonperturbative QCD. Currently, our experimental knowledge of the pion structure comes primarily from the Drell-Yan data for pion-nucleon/pion-nucleus scattering [1][2][3][4]. The valence quark distribution of the pion, q π v (x) with x being the fraction of the pion momentum carried by the active quark, has been extracted from these data [5][6][7]. Based on a next-to-leading order analysis including soft-gluon resummation, q π v (x) was found to behave as (1−x) 2 at large x [7]. On the other hand, theoretical predictions of q π v (x) have been made using various methods that are not fully consistent with this large-x behavior. For example, the parton model [8], perturbative QCD [9,10] and analysis from Dyson-Schwinger equations [11][12][13][14] suggest that the valence quark distribution should behave like (1 − x) a with a ≈ 2, while relativistic constituent quark models [15,16], Nambu-Jona-Lasinio models [17][18][19][20] and other arguments [21][22][23] favor a linear dependence on 1 − x at large x. For a review of the experimental and theoretical status of the pion parton distribution function (PDF), see Ref. [24]. Lattice QCD should be able to shed light on this puzzling disagreement, provided that its computational potential can be extended beyond the first few moments of PDFs.
This became possible recently due to the breakthrough made by large-momentum effective theory (LaMET) in direct lattice calculation of the x-dependence of PDFs [25,26]. According to LaMET, the full PDF, instead of its first few moments, can be directly ac-cessed from lattice QCD using the following method: 1) Construct an appropriate static-operator matrix element (now known as the quasi-PDF) that approaches the PDF in the large-momentum limit of the external hadron. The quasi-PDF constructed this way is usually hadronmomentum-dependent but time-independent, and, therefore, can be readily computed on the lattice. 2) Calculate the quasi-PDF on the lattice. 3) Convert it to the PDF through a factorization formula accurate up to power corrections suppressed by the hadron momentum. The existence of such a factorization is ensured by construction; for a proof, see Refs. [27][28][29].
In this paper, we carry out the first direct lattice calculation for the valence quark distribution of the pion using the LaMET approach. The calculation is done using clover valence fermions on an ensemble of gauge configurations with N f = 2+1+1 (degenerate up/down, strange and charm) flavors of highly improved staggered quarks (HISQ) [83] generated by the MILC Collaboration [84] with lattice spacing a = 0.12 fm, box size L ≈ 3 fm and pion mass m π ≈ 310 MeV. Our results are comparable quantitatively with the results extracted from experimental data [7] as well as from the Dyson-Schwinger equation [14].
From quasi-PDF to PDF in the pion: The quark PDF in the pion is defined as where the pion has momentum P µ = (P 0 , 0, 0, P z ), ψ f ,ψ f are the quark fields of flavor f , n µ = (1, 0, 0, −1)/ √ 2 is a lightlike vector, x denotes the fraction of pion momentum carried by the quark, and is the gauge link. The valence quark distribution is given by The quark quasi-PDF can be defined in a similar way to Eq. 1: except thatñ µ = (0, 0, 0, −1) is a spacelike vector with n · P = P z . As pointed out in Refs. [39,85], the Dirac matrix / n = γ z can also be replaced by γ t , which has the advantage of avoiding mixing with scalar PDF [47,57]. The choice of γ t will be used throughout this paper.
The factorization connecting the quasi-PDF and the PDF was first presented in Refs. [25,39] for bare quantities. Later on, it was shown [48,[54][55][56] that the renormalization factor of the quasi-PDF depends on the exponential of the Wilson line length times the Wilson line self energy counterterm. Nonperturbative renormalization can be carried out by extracting the counterterm from the heavy quark potential [37,54], using the RI/MOM scheme [29,35,46], or by forming a ratio of the lattice matrix elements of the quasi-PDF at two different momenta [29,76,[86][87][88]. Following our previous studies of nucleon PDFs [65], we perform nonperturbative renormalization in the RI/MOM scheme, but we also show the result from Wilson line renormalization in intermediate steps for comparison and to estimate the size of systematic uncertainties.
In the RI/MOM scheme, the bare coordinate-space matrix elementh(λñ) showing up on the right hand side of Eq. 3 can be renormalized nonperturbatively by demanding that the counterterm Z cancels all the loop contributions for the matrix element in an off-shell external quark state at a specific momentum [36,46]: and Then the nonperturbatively renormalized quasi-PDF can be matched to the PDF in the MS scheme: whereμ and µ denote the renormalization or cutoff scale for the quasi-PDF and the PDF, respectively. m π is the pion mass. The matching can be carried out perturbatively. At one loop, we define The matching kernel is the same as in Ref. [65]. Ideally, the continuum limit should be taken before matching such that lattice artifact can be removed and rotational symmetry can be recovered. However, only a single lattice spacing is used in this work.
For power corrections, the meson mass correction associated with the choice of Dirac matrix γ t is identical to that of the helicity distribution worked out in Ref. [31]. The O(Λ 2 QCD /(ñ · P ) 2 ) correction is numerically rather small in the present case. The renormalization and matching for the Wilson line renormalization scheme is shown in the Appendix.
Lattice calculation setup: In addition to the setup described in the Introduction, the gauge links are one step hypercubic(HYP)-smeared [89] with the clover parameters tuned to recover the lowest pion mass of the staggered quarks in the sea [90][91][92][93]. On these configurations, we calculate the time-independent, nonlocal (in space, chosen to be in the z direction) correlators of a pion with a finite-P z boost where U z is a discrete gauge link in the z direction, P = {0, 0, P z } is the momentum of the pion, and Γ = γ t .
To control the systematics due to contamination by excited states, we vary the Gaussian smearing parameter to best suppress the excited state, resulting in a clean ground-state pion. In addition, we use a simultaneous fit of the pion matrix element three-point correlators with four source-sink pion separations, 0.72, 0.84, 0.96 and 1.08 fm. First, we use the "two-state" analysis described in Ref. [93] to obtain the ground-state pion matrix elements using all 4 source-sink separations. A second extraction uses only the largest two separations. Finally, we use the "two-twoRR" analysis (see Ref. [93] for details), which includes an additional matrix element related to excited states. The detailed procedure for treating excited-state systematics can be found in Ref. [93] for the nucleon-charge case. Fig. 1 shows the real part of the matrix element for P z = 1.32 GeV using various combinations of data and analysis strategy. We will use matrix elements from the "two-twoRR" analysis in our PDF analysis. We use multiple values of pion momenta, P z = {0, 0, n 2π L }, with n ∈ {2, 3, 4}, which correspond to 0.86, 1.32 and 1.74 GeV, respectively.
Numerical results and discussions: Now we present our numerical results for the valence quark distribution in the pion and discuss their physical implications. We first Fourier transform the renormalized lattice data to momentum spacẽ form the valence distributionq R (x) +q R (−x), and then apply one-loop matching and meson mass corrections. The meson mass corrections are numerically rather small. To illustrate the impact of one-loop matching, we show the results before and after applying the matching at the largest momentum P z = 4 × 2π/L in Fig. 2. As can be seen from the plot, one-loop matching results in a sizable contribution and shifts of the original quasi-PDF towards the physical region [0, 1] in both schemes.
In our earlier work [32] on the nucleon PDF, we also proposed a "derivative" method to improve the truncation error in the Fourier transform in Eq. (9). To be concrete, we take the derivative of the renormalized nucleon matrix elements ∂ zhR (z), whose Fourier transform The pion boost momentum used here is Pz = 3×2π/L. The "two-two" and "two-twoRR" analysis uses all four source-sink separations while the "two-two2tsep" uses only the largest two source-sink separations. One can see that our final ground-state pion matrix elements are consistent among different ways of extracting the matrix elements. "two-twoRR" analysis includes additional excitedstate contribution; thus its errors are consistently larger even when comparing with those only extracting from larger separation where the three-point correlators are noisier.
where we can ignore boundary terms since the matrix elementh R (z) vanishes at large z. We also applied this method to the present lattice data. In Fig. 3, we show the results of the valence quark distribution for different pion momenta and renormalization schemes. For the RI/MOM result, we have chosen µ R = 3.7 GeV, p R z = 6 × 2π/L, and included statistical  as well as the systematic error of setting the unphysical scale p R z by varying it between 4 × 2π/L and 8 × 2π/L. For the Wilson line renormalization, only statistical errors are included since there is no extra unphysical scale in this scheme like p R z in RI/MOM. As can be seen from the figure, increasing P z tends to shift the distribution towards x = 0 and also lifts the peak at x = 0, but its impact is mild. Another important feature is that the RI/MOM result is consistent with 0 outside the physical region [0, 1] within errors, whereas the Wilson line renormalization one is not. This mainly reflects the importance of the higher-order matching kernel, since the one-loop matching in the two different schemes differ only by finite terms. We plan to derive higher-order matching, expecting that it will reduce the difference between the results in two different schemes.
It is worthwhile to stress that matching is a necessary step in converting the quasi-PDF to PDF. It yields sizeable contributions and changes, in particular for the distribution in the unphysical region. In Ref. [94], the authors studied the pion valence quasi-distribution using the Bethe-Salpeter wave function of the pion, and ob- GeV from RI/MOM scheme calculation (LP 3 ) , contrasted with analysis from Dyson-Schwinger equation [14] (DSE) at the scale 5.2 GeV and from a fit to Drell-Yan data in Ref. [7] (ASV) at 4 GeV.
served that for P z 2 GeV, by further increasing the pion momentum the quasi-PDF shrinks to the physical region very slowly. Actually we have observed a similar trend in our data. However, the matching plays an important role in reducing the contribution in the unphysical region, as can be seen from Fig. 2 above, but hasn't been taken into account in Ref. [94].
In Fig. 4, we compare our final result in RI/MOM scheme (LP 3 ) with computations from Dyson-Schwinger equation [14] (DSE) and from a phenomenological fit to Drell-Yan data [7] (ASV). We have set our renormalization scale to be µ = 4 GeV, in accordance with the experimental fit [7], whereas the DSE result is at 5.2 GeV. Outside the physical region, our result is consistent with 0. Within the physical region, our result decreases more slowly than the DSE and ASV results at large x, and has a lower peak around x = 0, as can be seen from the upper plot. This is expected to improve once we have lattice data at smaller pion masses. When plotted as xq π v (x), as was usually done in the literature, the discrepancy at small x gets suppressed, while it gets enhanced at large x.
We point out several potential sources of uncertainty or artifact in the above analysis, which we aim to improve in the future. First, the contribution at large x depends on the pion momentum as well as on the unphysical pion mass used in this calculation. If we have a larger pion momentum and a pion mass closer to its physical value, the contribution at large x will be further reduced, and accordingly, the small x contribution will be enhanced. Second, the matching implemented here is at one-loop order. It has a sizable effect and shifts the result towards the physical region. It is therefore important to investigate the impact of higher-order matching, in order to reduce uncertainties due to perturbative matching. This is also reflected by the difference between the results in two different schemes. Third, the present calculation is carried out at one lattice spacing, we'll need data at more lattice spacings to have a continuum extrapolation. Last but not least, we also need simulations with larger volumes to control the finite volume effect.

ACKNOWLEDGMENTS
We thank the MILC Collaboration for sharing the lattices used to perform this study. The LQCD calculations were performed using the Chroma software suite [95]. In this scheme, the nonperturbative renormalization readsh R (λñ) = Z 1 Z 2 e δmλh (λñ), The counterterm δm has been extracted nonperturbatively on the lattice with δm = 253(3) MeV at a = 0.12 fm [37,38]. The renormalization factors Z 1 and Z 2 come from the endpoint renormalization ofh(λñ) and can be fixed by an overall normalization. The matching kernel C (1) reads where Λ(x) = Λ 2 + x 2 P 2 z with Λ being a transverse momentum cutoff.