Valence parton distribution function of pion from fine lattice

TL;DR

This paper uses lattice QCD with Large Momentum Effective Theory to compute the valence parton distribution function of the pion, providing a first-principles determination with non-perturbative renormalization and matching.

Contribution

It introduces a lattice QCD approach with mixed action and RI-MOM scheme to extract the pion's valence PDF, testing the validity of 1-loop matching in this context.

Findings

Valence PDF of pion obtained from lattice QCD.

Validation of 1-loop matching procedure.

Consistent results across different lattice spacings.

Abstract

We present a lattice QCD study of valence parton distribution inside the pion within the framework of Large Momentum Effective Theory. We use a mixed action approach with 1-HYP smeared valence Wilson clover quarks on 2+1 flavor HISQ sea with the valence quark mass tuned to 300 MeV pion mass. We use $48^3 \times 64$ lattice at a fine lattice spacing $a=0.06$ fm for this computation. We renormalize the quasi-PDF matrix element in the non-perturbative RI-MOM scheme. As a byproduct, we test the validity of 1-loop matching procedure by comparing the RI-MOM renormalized quasi-PDF matrix element with off-shell quark external states as computed in the continuum 1-loop perturbation theory with the lattice results at $a=0.04$ and 0.06 fm. By applying the RI-MOM to ${\bar{\rm MS}}$ one-loop matching, implemented through a fit to phenomenologically motivated PDFs, we obtain the valence PDF of pion.

Equations

$$f(x)=\frac{1}{4\pi }\int d{\xi }^{-}{e}^{ix{P}^{+}{\xi }^{-}}⟨H(P)|\overline{\psi }({\xi }^{-}){\gamma }^{+}W({\xi }^{-},0)\psi (0)|H(P)⟩,$$

$$\tilde{q}(x,P_z)=\frac{1}{4\pi }\int dz{e}^{-ix{P}^{z}z}⟨H(P_z)|\overline{\psi }(z)\Gamma W(z,0)\psi (0)|H(P_z)⟩,$$

$$\tilde{q}(x,{\mu }_L,P_z)={\int }_{-1}^{+1}\frac{dy}{|y|}C\left(\frac{x}{y},\frac{y{P}_z}{\mu },\frac{{\mu }_L}{y{P}_z}\right)f(y,\mu ).$$

$$C_{2pt}^{s{s}^{\prime }}(t,P_z)=⟨\left[{\pi }_s^{+}(t,\vec{P}\right]{\left[{\pi }_{s^{\prime }}^{+}(0,\vec{P})\right]}^{†}⟩,$$

$$h_{\Gamma }^B(z,P_z)=\text{matrixelement}{E}_{\pi },P_z{\mathcal{O}}_{\Gamma }(z;\tau )E_{\pi },P_z,$$

$${\mathcal{O}}_{\Gamma }(z;\tau )=\sum _{\vec{x}}\left({\overline{u}}_x{W}_{x,x+\mathcal{L}}\Gamma u_{x+\mathcal{L}}-{\overline{d}}_x{W}_{x,x+\mathcal{L}}\Gamma d_{x+\mathcal{L}}\right),\text{where}$$

$$\Gamma ={\gamma }_t,{\gamma }_z,1;W_{x,x+\mathcal{L}}=\prod _{x^{\prime }=x}^{x+\mathcal{L}}{U}_3(x^{\prime }),$$

$$C_{\rm 3pt}^{SS}(t,\tau;z,P_{z})=\bigg{\langle}\big{[}\pi^{+}_{S}(t,\vec{P})\big{]}\mathcal{O}_{\Gamma}(z;\tau)\big{[}\pi^{+}_{S}(0,\vec{P})\big{]}^{\dagger}\bigg{\rangle},$$

$$R(t,\tau ;z,P_z,\Gamma )=\frac{C_{3pt}^{SS}(t,\tau ;z,P_z)}{C_{2pt}^{SS}(t,P_z)}.$$

$$R(t,\tau ;z,P_z,\Gamma )=\frac{{\sum }_{n,n^{\prime }}{A}_n{A}_{n^{\prime }}^{\ast }\text{matrixelement}{E}_n,P{\mathcal{O}}_{\Gamma }(z)E_{n^{\prime }},P)e^{-(E_{n^{\prime }}-E_n)\tau -E_{n}t}}{{\sum }_m|A_m{|}^2{e}^{-E_{m}t}},$$

$$⟨E_n,P_z|{\mathcal{O}}_{\Gamma }(z){|E_{n^{\prime }},P_z)⟩}^{\ast }={\varphi }_{\Gamma }⟨E_{n^{\prime }},-P_z|{\mathcal{O}}_{\Gamma }(z)|E_n,-P_z)⟩,$$

$$R_{\text{sum}}(t;z,\Gamma )=\sum _{\tau /a={\tau }_o}^{t/a-{\tau }_o}R(t,\tau ;z,P_z,\Gamma ).$$

$$R_{\text{sum}}(t;z,\Gamma )\simeq (t-2{\tau }_o)h_{\Gamma }^B(z,P_z)+\mathrm{const}+\mathcal{O}(e^{-(E_1-E_{\pi })t}).$$

$${\Lambda }_{\Gamma }(z,p)={⟨Q(p)⟩}^{-1}⟨\overline{u}(p){\mathcal{O}}_{\Gamma }(z,\tau )u(p)⟩{⟨Q(p)⟩}^{-1},$$

$$q_{\Gamma }(z,p)\equiv \mathrm{Tr}\left[\mathcal{P}{\Lambda }_{\Gamma }(z,p)\right],$$

$$q_{{\gamma }_t}^R(z,p,p^R)=Z_{{\gamma }_t{\gamma }_t}(z,p^R)q_{{\gamma }_t}(z,p),$$

$$q^{R}_{\gamma_{t}}(z,p,p^{R})\bigg{|}_{p=p^{R}}\equiv e^{ip^{R}_{z}z}.$$

$$q_{{\gamma }_z}^R(z,p,p^R)=Z_{{\gamma }_z{\gamma }_z}(z,p^R)q_{{\gamma }_z}(z,p)+Z_{{\gamma }_{z}1}(z,p^R)q_1(z,p).$$

$$\displaystyle q^{R}_{\gamma_{z}}(z,p,p^{R})\bigg{|}_{p=p^{R}}\equiv e^{ip^{R}_{z}z},$$

$$Z_{{\gamma }_z{\gamma }_z}(z,p^R)\mathrm{Tr}{\Lambda }_{{\gamma }_z}(z,p)+$$

$$\displaystyle\qquad Z_{\gamma_{z}1}(z,p^{R}){\rm Tr}\Lambda_{1}(z,p)\bigg{|}_{p=p^{R}}\equiv 0.$$

$$h_{{\gamma }_t}^R(z,P_z,p^R)$$

$$h_{{\gamma }_z}^R(z,P_z,p^R)$$

$$Z_q(p^R{)}^{-1}\frac{1}{12}\mathrm{Tr}\left(⟨Q(p^R){⟩}^{-1}{Q}_{\mathrm{tree}}(p^R)\right)=1,$$

$${\zeta }_{\Gamma }(z,p,p^R)=\frac{q_{\Gamma }^R(z,p,p^R)}{q_{\Gamma }^R(z,p,p)}-1,$$

$${\zeta }_{{\gamma }_t}(z,p,p^R)=\frac{Z_{{\gamma }_t{\gamma }_t}(z,p^R)-Z_{{\gamma }_t{\gamma }_t}(z,p)}{Z_{{\gamma }_t{\gamma }_t}(z,p)},$$

$$\displaystyle\zeta_{\Gamma}(z,p,p^{R})=\frac{\alpha_{s}C_{F}}{2\pi}\int_{-\infty}^{\infty}dx(e^{i(1-x)p_{z}}-1)\bigg{[}H\left(x,p\right)$$

$$\displaystyle-\left|\frac{p_{z}}{p_{R}}\right|H\left(1+\frac{p_{z}}{p_{z}^{R}}(x-1),p^{R}\right)\bigg{]},$$

$$q(z,p)=A{e}^{i\omega z}{e}^{-m_{\mathrm{scr}}|z|}{e}^{-c|z|},$$

$$m_{\mathrm{scr}}(z)-i\omega (z)\equiv -\frac{1}{a}\log \left(\frac{q(z+a,p)}{q(z,p)}\right)-c.$$

$$f_{u,d}(x)=Q_{u,d}(x)\theta (x)-Q_{\bar{u},\bar{d}}(-x)\theta (-x).$$

$$f_{u-d}(x)$$

$$f_v^{\pi }(x)$$

Full text

I Introduction

QCD factorization allows us to calculate the cross-section of hard hadronic processes in terms of the convolution of partonic cross-section and parton distribution functions Collins et al. (1989). Parton distribution for a hadron can be defined using hadronic matrix elements of appropriately chosen gauge invariant operators separated along the light cone. For example, the quark parton distribution function (PDF) of a hadron $H$ can be defined in terms of an operator bilocal in quark field $\psi$ as Soper (1977); Collins et al. (1989)

[TABLE]

where $W(\xi^{-},0)=\mathcal{P}e^{ig\int_{0}^{\xi^{-}}d\xi^{-}A^{+}}$ is the path-ordered straight Wilson Line on the light-cone, and the light-cone coordinates $\xi^{\pm}=(t\pm z)/\sqrt{2}$. A straight forward first principle calculation of PDF is not possible because lattice QCD is formulated in the Euclidean space-time, and thus, it cannot access quantities defined on the light-cone. To circumvent this problem, it has recently been proposed to calculate the quasi parton distribution function (qPDF), $\tilde{q}(x,P_{z})$, defined in terms of matrix elements of equal time, but spatially separated, quark bilinears Ji (2013) evaluated in a hadron state boosted to a large momentum $P_{z}$:

[TABLE]

where $\Gamma$ is either $\gamma_{z}$ or $\gamma_{t}$ for the unpolarized parton distribution addressed in this paper. Here, $W(z,0)$ is a straight spatial Wilson line joining the quark and anti-quark. For sufficiently boosted hadrons, one can use the Large Momentum Effective Theory (LaMET) Ji (2014) to relate the qPDF to PDF through a convolution with a matching kernel $C$ as

[TABLE]

Here $\mu_{L}$ and $\mu$ are the renormalization scales of the schemes in which the qPDF and PDF are defined. For the latter, ${\overline{\rm MS}}$ scheme is used and $\mu$ is referred to as the factorization scale. The matching kernel is perturbative and hence universal for all the hadrons. Therefore, it is calculated using quark external states in a chosen gauge. Such calculations at 1-loop order have been performed using the cutoff scheme Xiong et al. (2014) as well as in the ${\overline{\rm MS}}$ scheme Constantinou and Panagopoulos (2017); Liu et al. (2018); Stewart and Zhao (2018). There are also related approaches to calculate the PDF from the lattice that use similar logic but differ in details, like the pseudo-PDF approach proposed in Ref. Radyushkin (2017); Orginos et al. (2017) and the use of good lattice cross sections Ma and Qiu (2018a, b). The latter includes the current-current correlators Sufian et al. (2019).

Using LaMET and related approaches, various attempts have been made to calculate the unpolarized and polarized iso-vector quark distribution of the nucleon Liu et al. (2018); Lin et al. (2018a); Chen et al. (2018a); Alexandrou et al. (2018a, b). The first studies of the valence quark distribution for the pion have also been presented Sufian et al. (2019); Chen et al. (2018b). One important issue in the calculation of the PDF from the lattice is the renormalization and matching. As indicated above, the PDF and qPDF are usually defined in different renormalization schemes. The qPDF, which is calculated on the lattice, needs a non-perturbative renormalization scheme because of the self-energy divergence of the Wilson line Ishikawa et al. (2016), and this is usually implemented using the RI-MOM scheme Chen et al. (2018c) defined using external off-shell quark states accessible on the lattice. Then, one has to match the qPDF in this lattice renormalization scheme to the PDF in the ${\overline{\rm MS}}$ scheme through Eq. (3). This is achieved through the convolution using the matching kernel between the RI-MOM and ${\overline{\rm MS}}$ schemes that is perturbatively calculated in the continuum theory using dimensional regularization Stewart and Zhao (2018). One could also define the qPDF operator in the ${\overline{\rm MS}}$ scheme and then perform the matching between PDF and qPDF Constantinou and Panagopoulos (2017); Alexandrou et al. (2017). The current status of this field, including the comparison with the phenomenological PDF and the issue of renormalization, is reviewed in Refs. Zhao (2019); Cichy and Constantinou (2019); Monahan (2018).

In principle, Eq. (3) offers a way to calculate PDF from the lattice, but it is unclear as to what extent this is actually feasible given the various assumptions that go along with the equation implicitly. For example, at any finite hadron momentum $P_{z}$, the Eq. (3) suffers from ${\cal O}\left(\Lambda^{2}_{\rm QCD}/(x^{2}P^{2}_{z})\right)$ higher twist corrections. This is closely related to the assumption that the perturbative calculation, currently truncated at 1-loop order, is able to capture the renormalization as well as the matching of the qPDF matrix element over a range of quark-antiquark separations, $z$ – to be in the perturbative regime, one would expect $z$ to be smaller than or about $\mathcal{O}(1)$ fm. It is also important to ensure that $aP_{z}<1$ to make sure we are not overcome with lattice artifacts Xiong et al. (2017). Therefore, a closer look at this new methodology is warranted and is actively being studied Alexandrou et al. (2019); Liu et al. (2018). The aim of this paper is to explore these issues further by using finer lattices than what is being used in the qPDF literature, and use pion as a case study. The smaller mass of the pion makes it easier to achieve a large boost, the numerical calculations are expected to be less expensive and it also helps suppress the target mass correction by ensuring $m_{\pi}\ll P_{z}$. We focus on the valence PDF of pion since it can be accessed using the isotriplet $u-d$ PDF, and thereby, avoid mixing with the gluon sector. In our study, we will use the renormalization and matching strategy outlined in Chen et al. (2018c); Liu et al. (2018). The pion valence PDF has been determined through a leading order and next to leading order analyses of the experimental data Badier et al. (1983); Betev et al. (1985); Conway et al. (1989); Owens (1984); Sutton et al. (1992); Gluck et al. (1992, 1999); Wijesooriya et al. (2005); Aicher et al. (2010), but it is much less constrained than the nucleon PDF and therefore, the lattice calculations may have more impact in this case, especially in constraining the $x\to 1$ limit which is not yet well established.

The paper is organized as follows. In section II, we discuss our lattices setup. In section III, we present the calculations of the two point function of the boosted pion and check how reliable the extractions of the ground state and the first excited state are. In Section IV, we present our results for the pion three point function that defines the qPDF. Here, we also discuss the problem of excited state contamination. In section V, we discuss the non-perturbative renormalization as well as the validity of 1-loop matching. Our results on the renormalized pion qPDF and the matching to PDF are presented in section VI. Some technical aspects of the calculations are discussed in the Appendices. Preliminary results on this work have been reported in conference proceedings Petreczky et al. (2018); Karthik et al. (2018); Shugert et al. (2018).

II Lattice setup

We performed the calculations of the pion two-point and three-point functions needed to obtain the qPDF using the Wilson-Clover action for valence quarks on 1-HYP smeared gauge configurations Hasenfratz and Knechtli (2001) and the Highly Improved Staggered Quark (HISQ) action Follana et al. (2007) in the sea. We used the 2+1 flavor gauge configurations corresponding to lattice size $48^{3}\times 64$ and the lattice spacing of $a=0.06$ fm generated by the HotQCD collaboration Bazavov et al. (2014). In addition to this ensemble, we also used $64^{4}$ HISQ lattices Bazavov et al. (2014) with the lattice spacing $a=0.04$ fm for the study of the non-perturbative renormalization (NPR). In both the ensembles, the sea quark mass was tuned to a pion mass of 160 MeV. A similar setup was used by the PNDME collaboration albeit for 2+1+1 flavor MILC configurations (c.f., Ref. Liu et al. (2018)). For the valence quark masses, we used the values $am=-0.0388$ (i.e., $\kappa=0.12623$) for the $a=0.06$ fm ensemble and $am=-0.033$ (i.e., $\kappa=0.12604$) for the $a=0.04$ fm ensemble, which are tuned such that the pion mass, $m_{\pi}$, is 300 MeV. We did not see any exceptional configurations for these valence quark masses in our calculations.

We used higher statistics at smaller quark-antiquark separations $z$ than at larger ones; to be exact, we used $216$, $100$ and $48$ gauge configurations for $|z|/a\in[0,8]$, $(8,16]$ and $(16,24]$ respectively. We further improved the statistics by using the All-Mode Averaging (AMA) Shintani et al. (2015) technique in the computations of the two- and three-point functions, with 32 sloppy calculations to one exact solve for each configuration. For the exact and sloppy inversions, we used the stopping criterion of $10^{-10}$ and $10^{-4}$, respectively. In our study, we will consider the valence quark distribution, which in turn is related to the iso-vector $u-d$ quark distribution in the pion, and thus we do not compute the quark line disconnected diagrams.

For a reliable extraction of qPDF, a good overlap of the source operator with the pion state is necessary so as to project out the ground state at as small source-sink separation as possible. The quark sources with Gaussian profile, typically implemented through a gauge covariant Wuppertal smearing Gusken et al. (1989), are used for this purpose when the hadron is at rest. However, for the fast moving hadrons that are required in the qPDF framework, the use of the Gaussian sources is no longer sufficient and this necessitates the usage of the boosted Gaussian sources Bali et al. (2016) instead. Since we are interested in the calculation of the pion two and three point function at several values of the pion momenta and several source-sink separations, we found it more practical to implement the Gaussian sources by using the Coulomb gauge instead of implementing the Wuppertal smearing. We found the optimal size of the Gaussian profile to be about 0.3 fm, which roughly corresponds to 90 steps of Wuppertal smearing. We checked that in terms of the signal-to-noise ratio, the Wuppertal and Coulomb-gauge Gaussian sources are similar (see Appendix E). In the next section, we discuss the boosted sources in detail and the energy levels of the boosted pion. In Appendix A, we have explained the construction of boosted sources in detail.

Out of 216 gauge configurations used in our calculations, 24 gauge configurations were analyzed using the GPU cluster in BNL to calculate two-point and three-point correlation functions. These calculations were performed entirely on GPU using the QUDA suite Clark et al. (2010); Babich et al. (2011); Clark et al. (2016), including the inversion of the fermion operator with multigrid algorithm, communication between GPU devices to perform covariant shifts, and the necessary spin-color matrix multiplications. In QUDA, the communications between GPUs on the same node are implemented through MPI or as direct peer-to-peer communications between the GPU devices. We have found that on rare occasion the QUDA peer-to-peer communications did not finish by the time the computations started. These rare glitches happened randomly. We checked, however, that these glitches did not affect our results noticeably compared to other errors.

III Two point function of the boosted pion

We calculated the two point functions of the positively charged pion ($\pi^{+}=\overline{d}u$),

[TABLE]

for a spatial pion momentum $\vec{P}=(0,0,P_{z})$ which is non-zero only along the $z$-direction, using the pion source and sink $\pi^{+}_{s}(0,\vec{P})$ and $\pi^{+}_{s^{\prime}}(t,\vec{P})$, respectively. The values of momenta in lattice units are $aP_{z}=\pm 2\pi n_{z}/48$ for $n_{z}$ ranging from 0 to 5, which in physical units correspond to $P_{z}=0,0.43,0.86,1.29,1.72$ and $2.15$ GeV, respectively. We always used the Coulomb gauge Gaussian smeared-source ($s=S$), and either a smeared sink ($s^{\prime}=S$) or point sink ($s^{\prime}=P$). In the rest of the paper, we will refer to the smeared-source and smeared-sink set-up to be SS, and we will refer to the smeared-source point-sink set-up as SP.

For the lowest two momenta, we used the usual Gaussian sources. To improve the signal for the higher momenta, we followed Ref. Bali et al. (2016) and used boosted sources in which the valence quarks are boosted to a momentum $k_{z}=\zeta P_{z}$, with $\zeta$ being a tunable parameter. Naively, one might expect that the optimal choice to be $\zeta=0.5$. However, we found that the optimal choice of $\zeta$ for the pion in terms of the signal-to-noise ratio is between $0.6-0.75$. For $P_{z}=0.86$ GeV the signal-to-noise ratio is not very sensitive to the value of $\zeta$. These findings are in agreement with Ref. Bali et al. (2016). We discuss the optimization of boosted sources further in Appendix E. Since, we need to create a source for each value of $\zeta$, we used $k_{z}=2(2\pi/48)$ for $n_{z}=2,3$ and $k_{z}=3(2\pi/48)$ for $n_{z}=4,5$, corresponding to the choices of the parameter $\zeta=1,2/3,3/4$ and $3/5$ for $n_{z}=2,3,4$ and $5$ respectively. We have shown the corresponding effective masses for the SS two point functions in Fig. 1. By using the boosted smeared sources, one can see that a reasonable signal for the two point correlation function can be obtained up to source-sink separations $t=12a$ for all momenta except for the highest momentum $P_{z}=2.15$ GeV. Simply from the data points in Fig. 1, we see that the effective masses approach a plateau corresponding to the continuum dispersion relation $E_{\pi}(P_{z})=\sqrt{P_{z}^{2}+m_{\pi}^{2}}$, shown as the horizontal lines. The effective mass approaches the plateau region at larger source-sink separations when the momentum is increased, as one would expect from the shrinking gap between the ground and excited states as the pion is boosted.

As we will discuss next, we used source-sink separations $t=8a,10a$ and $12a$ for the computation of three-point functions. Therefore, we needed to analyze the excited state contribution to the SS and SP correlators to perform the infinite source-sink extrapolations. For this, we performed multi-exponential fits on the SS and SP pion two point functions in the interval $t\in[t_{\rm min},t_{\rm max}]$ in order to extract the energy levels. For fixed $t_{\rm min}$, we varied $t_{\rm max}$ and checked the sensitivity of the result to $t_{\rm max}$. Then, we repeated the procedure for different values of $t_{\rm min}$. We found that we were able to reliably extract the ground state $E_{\pi}(P_{z})$ as well as the first excited state $E_{1}(P_{z})$ using the four-parameter two-state fits to the SP correlator instead of using the SS correlator. This could be due to the fact that the contributions from the high-lying energy levels are smaller in the SP correlator compared to the SS correlator stemming from the possible cancellations between the positive as well as the negative amplitudes that are allowed in the SP correlator. In the top-left and bottom-left panels of Fig. 2, we have shown the systematics of the two-state fits to the SP correlator. In the top-left and the bottom-left figures, we have shown the dependence of the best fit values of $E_{\pi}$ (blue circles) and $E_{1}$ (black circles) as a function of $t_{\rm min}$ used in the fits. For a given $t_{\rm min}$, the data points from two values of $t_{\rm max}$ have been clubbed together for $P_{z}=0$ and 1.29 GeV respectively, and it demonstrates that there is no dependence on $t_{\rm max}$. The ground state is seen to compare well with the expectation from the dispersion relation shown by the black solid lines. The red band shows the values of $E_{1}$ chosen as the best estimate of the first excited state.

On the top-right and bottom-right panels of Fig. 2, we show similar plots for the first excited state $E_{1}$ as estimated using the SS correlator. The statistical errors of the excited state energy $E_{1}$ in the simple two state fits (magenta diamond) quickly grow large with increasing $t_{\rm min}$ and thus, these fits turned out to be of limited use. Therefore, we performed constrained two- and three-state exponential fits for the SS correlator with the ground state energy fixed to $E_{\pi}=\sqrt{P_{z}^{2}+m_{\pi}^{2}}$ with $m_{\pi}=300$ MeV, and imposing a prior on $E_{1}$ using its best estimate from the SP correlator – that is, we added the term $(E_{1}-E_{1,{\rm prior}})^{2}/\sigma_{\rm prior}^{2}$ to the $\chi^{2}$ with $E_{1,{\rm prior}}$ and $\sigma_{\rm prior}$ being the mean and error of $E_{1}$, respectively, as determined from the SP correlator. The $t_{\rm min}$ dependence of the resulting $E_{1}$ from the constrained two-state fit (black circles) and constrained three-state fit (blue triangles) are shown in the top-right and bottom-right panels. The two-state fits of the SS two-point correlator largely overestimate the energy of the first excited state for small $t_{\rm min}$ whether or not priors are used, and there is a significant dependence on $t_{\rm min}$. One should use $t_{\rm min}\geq 6a$ to obtain reliable results for the first excited state from the SS correlator. The three-state fits with priors on $E_{\pi}$ and $E_{1}$ give energies of excited states that are the same within errors for the SS two point correlators and show almost no $t_{\rm min}$ dependence. In summary, we determined the lowest three energy levels using the SP correlator and then determined the corresponding amplitudes $|A_{n}|=|\bra{0}\pi^{+}_{S}(0,\vec{P})\ket{E_{n},P_{z}}|$ of these excited states in the SS correlator through a constrained fit analysis.

In Fig. 3, we show the three energy levels obtained from the different fits discussed above, as a function of $P_{z}$. For $P_{z}=0$, we compare our result with the energy levels that would correspond to the pion resonances $\pi(1300)$ and $\pi(1800)$ from the PDG Tanabashi et al. (2018). In order to account for the 300 MeV pion mass, we shifted the PDG values by 0.161 GeV as an approximation and these are shown as the two arrows in Fig. 3. Our estimate of the first excited state energy agrees with this shifted mass of $\pi(1300)$. We also show the expected $P_{z}$-dependence of $E_{1}(P_{z})$ assuming a particle-like dispersion and this describes the actual data very well. The energy of the second excited state is much larger than expected, meaning that the third state effectively parametrizes several higher lying states. As one can see from the figure the energy gap between $E_{\pi}$ and $E_{1}$ shrinks with increasing $P_{z}$ as expected.

The results on the excited state energies will be important for the analysis of the pion three-point function discussed in the next section.

IV Extraction of the bare quasi-PDF matrix elements from the three-point functions

The next step is the calculation of the bare qPDF matrix element

[TABLE]

where the bilocal $u-d$ qPDF operator in a time-slice $\tau$ involving a quark and an antiquark separated along the $z$-direction by ${\cal L}=(0,0,0,z)$ is given by

[TABLE]

and it is made gauge-invariant by the Wilson line $W_{x,x+{\cal L}}$. The Dirac $\gamma$ matrices in the qPDF operator are in the Minkowskian convention. The state $\ket{E_{\pi},P_{z}}$ denotes the onshell ground state pion with momentum $P_{z}$. In addition to the natural choices of $\Gamma=\gamma_{t}$ and $\gamma_{z}$ that approach $\gamma^{+}$ in the light-cone limit, we also considered $\Gamma=1$. This choice of $\Gamma$ is needed because under renormalization, $\mathcal{O}_{\gamma_{z}}(z)$ mixes with $\mathcal{O}_{1}(z)$ Constantinou and Panagopoulos (2017). We applied one-level of HYP smearing to the links entering $W_{x,x+{\cal L}}$ in order to reduce the noise. Since the qPDF calculation involves values of $z\sim{\cal O}(a)$, we checked that there is no significant difference between the renormalized matrix elements using the smeared and unsmeared Wilson line. To obtain the bare matrix element $h^{B}_{\Gamma}(P,z)$, we computed the three-point function at different source-sink separations $t$ and operator insertion point $\tau$

[TABLE]

and constructed the ratio of the three-point function to two-point function,

[TABLE]

The reader can refer to Appendix B for a detailed description of the construction of three-point functions. The two-point function is always real when the source and sink are of the same type. The three-point function for the $u-d$ qPDF operator ${\cal O}_{\Gamma}(z)$ in a pion external state is real at all $z$ for $\Gamma=\gamma_{z},\gamma_{t}$ and purely imaginary for $\Gamma=1$ (refer Appendix C). Inserting a complete set of states in the above equation,

[TABLE]

with $E_{n+1}\geq E_{n}$, and $E_{0}=E_{\pi}$. It is easy to see that in the infinite $t$ limit, $R(t,\tau;z,P_{z},\Gamma)$ is equal to $h^{B}_{\Gamma}(z,P)$. The above equation holds for the infinite time extent. For a finite time extent, the effects of periodic boundary condition should be taken into account. This turns out to be important for $P_{z}=0$, while for non-zero $P_{z}$ the effect is negligible as discussed in Appendix G. In practice, one truncates the sums in Eq. (10) at some value $n$, and then obtain $h^{B}_{\Gamma}(z,P_{z})$ by fitting the $t$ and the $\tau$ dependence of $R(t,\tau;z,P_{z},\Gamma)$ using $\bra{E_{n},P}\mathcal{O}_{\Gamma}(z)\ket{E_{n^{\prime}},P)}$ as fit parameters. In the fits, the values of $A_{n}$ and $E_{n}$ were held fixed at values determined from the two-state fit analysis on the SS correlators. In what follows, we will refer to this method of fitting using $n$-state Ansatz to the data between $\tau/a>\tau_{o}$ and $\tau/a<t/a-\tau_{o}$ as Fit($n$,$\tau_{o}$). In this method, the excited states are suppressed by $\exp(-(E_{n}-E_{\pi})t/2)$. For $z=0$, it is easy to see that $R(t,\tau;z,P_{z},\Gamma)$ is symmetric in $\tau$ around the mid-point $\tau-t/2$. For $z\neq 0$ and $P_{z}\neq 0$, we only have the following relation (see Appendix D),

[TABLE]

where $\phi_{\Gamma}=1$ for $\Gamma=\gamma_{t},1$ and $\Phi_{\Gamma}=-1$ for $\Gamma=\gamma_{z}$. Thus, generically the matrix elements $\bra{E_{n},P_{z}}\mathcal{O}_{\Gamma}(z)\ket{E_{n^{\prime}},P_{z})}$ and $\bra{E_{n^{\prime}},P_{z}}\mathcal{O}_{\Gamma}(z)\ket{E_{n},P_{z})}$ are independent and the number of fit parameters is thereby increased. Based on the above relation, we constructed appropriate averages using both the positive and negative values of momenta to increase the statistics. However, in practice the gain was marginal.

We demonstrate the extraction of the matrix element using the fit method in Fig. 4, where the ratios $R(t,\tau;z,P_{z},\Gamma)$ are shown for $\Gamma=\gamma_{t}$ and $\gamma_{z}$ qPDFs for the $P_{z}=1.29$ GeV pion. The results on $R(t,\tau;z,P_{z},\Gamma)$ and the extraction of matrix elements for the other values of $P_{z}$ are given in Appendix G. In the figure, we show the data at $t/a=8,10$ and $12$ along with the result of Fit(2,2). Using the fit, the results for the $t\to\infty$ extrapolations are shown with the horizontal bands. We also performed the three states fit of $R(t,\tau;z,P_{z},\Gamma)$ and the picture looks similar. In this case, the data points at all $\tau-t/2$ could be described by the fit. The $t\to\infty$ extrapolations from the three state fit gave results consistent with the two state ones, albeit with larger errors. A closer look at Fig. 4 (and also from Fig. 21) reveals that the excited state contribution is larger for $\Gamma=\gamma_{z}$ than for $\Gamma=\gamma_{t}$. Furthermore, the excited state contribution grows with increasing $z$. The non-symmetric nature of $R(t,\tau;z,P_{z},\Gamma)$ for $z\neq 0$ is also apparent in the figure. We expect that $h_{\Gamma}(z=0,P)=1$ for $\Gamma=\gamma_{t}$ because of the charge conservation once a proper renormalization is implemented and the continuum limit is taken. Our extrapolation procedure gives a result for the bare matrix element which is larger than one at all the values of momenta as can be seen in Fig. 4, as well as from Fig. 22 in Appendix G where, in addition, one can also see that $h^{B}_{\gamma_{t}}(z=0,P_{z})$ is independent of $P_{z}$. Thus, any deviation of $h^{B}_{\gamma_{t}}(z=0,P_{z})$ away from unity should be taken care of by the renormalization. We will see in Section V that this is indeed the case.

Alternatively we can use the summation method Maiani et al. (1987) to obtain $h^{B}_{\Gamma}$. Here one sums over all $\tau/a$ minus a certain number of end points $\tau_{o}$

[TABLE]

We will refer to this method as sum($\tau_{o}$). For large $t$, one would find a linear behavior in $t$ of $R_{\rm sum}$ as

[TABLE]

The advantage of this method is that the excited state contributions are suppressed as $\exp[-(E_{n}-E_{\pi})t]$ instead of being suppressed as $\exp[-(E_{n}-E_{\pi})t/2]$ in the fitting method. We show a sample result using sum(1) and sum(2) in Fig. 5 for $\Gamma=\gamma_{z}$ and $z=0$. We see that $R_{\text{sum}}(t;z,\gamma_{z})$ can be well fitted by a straight line in $t$, and the slope gives the value of the matrix element. As a cross-check, we also show the expected curve for $R_{\text{sum}}(t;z,\Gamma)$ using our best fit from Fit(2,2) as the dashed curves. It can be seen that the difference between a simple straight line fit and the curve from Fit(2,2) is small. One can also note that sum(1) and sum(2) are almost parallel, meaning that the extracted matrix element is independent of $\tau_{o}$ confirming that the method works well.

To better understand the systematic effects due to excited state contaminations, one can look at the case $z=0$ in detail, where the statistical errors are the smallest. The bare matrix element $h^{B}_{\gamma_{z}}(z=0,P_{z})$ after renormalization is expected to be proportional to the hadron velocity, $P_{z}/E_{\pi}(P_{z})$. One can take the ratio of matrix elements $h^{B}_{\gamma_{z}}(z=0,P_{z})/h^{B}_{\gamma_{t}}(z=0,P_{z})$ to avoid the issues of renormalization. The results for this ratio of matrix elements is shown in Fig. 6 along with the curve for $P_{z}/\sqrt{P_{z}^{2}+m_{\pi}^{2}}$. We see that our lattice results follow the expectations reasonably with small $3-4\%$ deviations from the expected result at small $P_{z}$. A reason for this could be the large systematic uncertainty in $\gamma_{z}$ matrix element due to the relatively larger excited state extrapolations required. We see that within errors, the two-state fit, three-state fit and the summation methods are consistent.

In Fig. 7, we show the results for $h^{B}_{\Gamma}(z,P_{z})$ as functions of $z$ for the two highest momenta $P_{z}=1.29$ and $1.72$ GeV determined using the HYP smeared Wilson line. Since the real part is symmetric about $z=0$, we have only shown the data for $z\geq 0$. At each $z/a$, we have shown the resulting $t\to\infty$ extrapolated results using Fit(3,1), Fit(2,2), sum(2) and sum(2) methods, and these points at a given $z/a$ are slightly displaced for better visibility. We see that the results from all these methods agree with each other within the errors. For $\gamma_{z}$, some tension between the summation method, and the two and three state fits is observed at larger $|z|$. At larger values of $z$ the matrix elements are suppressed partly because of the larger value of $P_{z}$ and in part by the divergent self-energy contribution in the spatial Wilson line. The latter will be removed upon renormalization as we will see in the next section. Having demonstrated a robust determination of the matrix element using multiple $t\to\infty$ strategies, we will use the matrix elements obtained using Fit(2,2) in the rest of the paper.

So far we discussed results on the three-point function obtained using 1-HYP smearing for the spatial link. We also performed calculations using the unsmeared spatial Wilson line. In this case, the bare matrix element rapidly decreased with $z$ due to the larger value of the Wilson line self-energy divergence. However, we found that the results, after non-perturbative renormalization (discussed next in Section V), were similar to those obtained with smeared Wilson lines within errors. The main difference between the renormalized three-point function obtained with the smeared and unsmeared Wilson line is that, for the latter the statistical errors at large $z$ are significantly larger.

V RI-MOM non-perturbative renormalization and its comparison with 1-loop

In the last section, we discussed the extraction of the bare qPDF matrix element which has to be renormalized. The renormalizability of qPDF has been recently demonstrated to all orders of perturbation theory Ishikawa et al. (2017); Ji and Zhang (2015). In addition to the quark wavefunction renormalization $Z_{q}$ and the composite operator renormalization required for $z=0$, the qPDF operator at non-zero $z$ requires additional renormalization due to the the UV divergence present in the Wilson line connecting the quark and anti-quark Polyakov (1980). When a lattice fermion that breaks chiral symmetry at finite lattice spacings is used, as is the case in this paper, it has been shown that only the renormalization of $\gamma_{t}$ qPDF operator is purely multiplicative, while the $\gamma_{z}$ qPDF operator mixes with the scalar qPDF Chen et al. (2018c); Constantinou and Panagopoulos (2017). A renormalization scheme that is implementable on the Euclidean lattice is the RI-MOM scheme and it is now standard in the lattice QCD literature. The corresponding RI-MOM counter-term for the qPDF operator in leading order perturbation theory has been worked out using off-shell quark external states Stewart and Zhao (2018), and it is one of the ingredient used in the perturbative matching of the RI-MOM renormalized qPDF to the ${\overline{\rm MS}}$ PDF. In this section, we discuss the renormalization procedure, and then compare the running of the renormalization constants as determined on the lattice with the corresponding perturbative expectations. This allows us to quantitate the validity of the leading order perturbation theory and matching.

For non-perturbative renormalization, we compute the expectation value of qPDF operator between offshell quark external states with momentum $p$. We refer to the momentum of quark in the direction of Wilson-line as $p_{z}$ and the magnitude of the component perpendicular to the Wilson-line as $p_{\perp}$. For these computations, we use Landau gauge fixing. Let $\Lambda_{\Gamma}(z,p)$ be the quark-line amputated bare qPDF,

[TABLE]

where $Q(p)$ is the quark propagator $\langle u(p)\overline{u}(p)\rangle$ and $u(p)=\sum_{x}u_{x}e^{-ip.x}$. Let us define the bare qPDF after projection as

[TABLE]

consistent with the definition used in perturbative calculations. Here, ${\cal P}$ is the operator used to project onto one of the $\gamma$-matrices $\Gamma=\gamma_{\alpha}$, and ${\rm Tr}(\ldots)$ is a trace over both color and Dirac indices. Based on previous works Stewart and Zhao (2018); Chen et al. (2018c), we will use $\not{p}$-projection for which ${\cal P}=\not{p}/\left(12p_{\alpha}\right)$. Alternatively, one can use ${\mathcal{P}}=\Gamma$ Alexandrou et al. (2017) or the minimal projection Stewart and Zhao (2018). In the case of $\Gamma=\gamma_{t}$, since the renormalization is simply multiplicative, the renormalized quark qPDF is given by

[TABLE]

where the $z$-dependent RI-MOM renormalization constant $Z$ is determined using the renormalization condition set at momentum $p^{R}$ as

[TABLE]

The right-hand-side of the above equation is the tree-level value of $q_{\gamma_{t}}$. The renormalization constant so obtained is in general a complex number. For $\Gamma=\gamma_{z}$, we have to take care of mixing with the scalar $\Gamma=1$. Hence, the renormalized qPDF is defined as

[TABLE]

The diagonal part $Z_{\gamma_{z}\gamma_{z}}$ and the mixing term $Z_{\gamma_{z}1}$ are determined using the two RI-MOM conditions Chen et al. (2018c)

[TABLE]

Using the renormalization constants $Z$ determined above using quark external states, the renormalized pion qPDF can also be determined by

[TABLE]

where $Z_{q}$ is the quark renormalization, that can be determined using the condition Alexandrou et al. (2011)

[TABLE]

where $Q(p)$ is the quark propagator determined using the Landau gauge and $Q_{\rm tree}$ is the free quark propagator for which we use the free massless Wilson-Dirac propagator.

In Fig. 8, we show the renormalization factors using the above RI-MOM renormalization conditions on the $0.06$ fm ensemble. On the top left panel of Fig. 8, we show the real and imaginary parts of $Z_{\gamma_{t}\gamma_{t}}$ determined at $p_{z}^{R}=1.29$ GeV and $p_{\perp}^{R}=1.49$ GeV. The rapid, almost exponential, increase in $Z$ with $z$ is due to the self-energy divergence present in the bare Wilson line that connects the quark and antiquark in the qPDF operator. This divergent piece, $e^{c|z|}$, cannot be captured perturbatively and it needs to be determined nonperturbatively in a particular scheme. However, this might not be an issue for the one-loop matching if $e^{c|z|}$ cancels exactly between the renormalization factors and the bare qPDF operator. Therefore, we remove $e^{c|z|}$ from the renormalization constant that is shown in the top-left panel, and display the result in the top-right panel. The value of $c$ for our $a=0.04$ fm ensemble was determined in Bazavov et al. (2018), and for 1-HYP Wilson line $ca=0.1586$. This removal of Wilson line self-energy reduces the almost exponential dependence of $Z(z)$ to a weak dependence on $z$. In fact, we see that both the real and imaginary parts of $Z_{\gamma_{t}\gamma_{t}}$ remain ${\cal O}(1)$ even up to $z=1$ fm, thereby providing a qualitative justification for the usage of leading order perturbation theory to describe the lattice data at short $z$ and at high quark momenta. We show similar data in the bottom left and right panels for the Z-factors for $\gamma_{z}$ qPDF. In this case, we have the diagonal factor $Z_{\gamma_{z}\gamma_{z}}$ as well as the off-diagonal factor $Z_{\gamma_{z}1}$ to take care of mixing with scalar on the lattice. We show $Z_{\gamma_{z}\gamma_{z}}$ and $Z_{\gamma_{z}1}$ as the filled and unfilled symbols in the bottom panels respectively. We observe that the imaginary part of $Z_{\gamma_{z}\gamma_{z}}$ is small compared to the real part. This is not the case for $Z_{\gamma_{t}\gamma_{t}}$, which in turn will affect the asymmetry of the $u-d$ qPDF $\tilde{q}_{u-d}(x)$ of pion about $x=0$. We also note that the mixing of $\gamma_{z}$ with the scalar is a minor 5-10% effect, but we nevertheless take care of it in our calculation.

As we discussed in the last section, the matrix element at $z=0$, $h_{\gamma_{t}}$ is the local current operator which will be exactly conserved in the continuum limit. Hence, $Z_{q}Z_{\gamma_{t}\gamma_{t}}(z=0)$ is the vector current renormalization factor $Z_{V}$ and the dependence of $Z_{V}$ on $p$ will give us an idea of the leading $(pa)^{2}$ perturbative lattice artifacts for values of $p\gg\Lambda_{\rm QCD}$ as well as the other higher order (or perhaps non-perturbative) contributions to this lattice correction to $Z_{V}$ at smaller renormalization scales Alexandrou et al. (2011). In Fig. 9, we show $Z_{q}$ determined using Eq. (25), the value of $Z_{\gamma_{t}\gamma_{t}}$ at $z=0$ as well as their product $Z_{V}$ as a function of $(pa)^{2}$. One sees a reasonable plateau for $Z_{V}\approx 0.97$ only for $(pa)^{2}>2$. For comparison, the value of $Z_{V}$ as obtained from the bare pion isospin charge $h^{B}_{\gamma_{t}}(z=0,P_{z}=0)$ is 0.961(3). The values of $Z_{V}$ determined from $h^{B}_{\gamma_{t}}(z=0,P_{z})$ at the other non-zero $P_{z}$ also give consistent values. With the uncertainties of choosing the scaling region in $(pa)^{2}$ to take the $(pa)\to 0$ limit of $Z_{V}$, we expect the $Z_{V}$ to be in the range 0.97 to 0.99. For relatively smaller values of renormalization momenta $(pa)\approx 1-1.5$, chosen such that the renormalization scales lie in the vicinity of the pion momenta used in this paper, one sees noticeable, but small $~{}5\%$ dependence on $pa$. We used the value of $Z_{q}$ estimated at the same value of $p$ as used in $Z_{\gamma_{t}\gamma_{t}}$ for renormalizing our pion qPDF.

V.1 Comparison with leading order perturbation theory for $z<0.3$ fm

We will now investigate in a quantitative way the agreement/disagreement of the lattice determination of the RI-MOM renormalized amputated quark qPDF at $z<0.3$ fm which one can expect to be in the perturbative regime. For this, we construct a quantity $\zeta_{\Gamma}(z,p,p^{R})$ in the following way

[TABLE]

where $q^{R}_{\Gamma}(z,p,p)=e^{ip_{z}z}$ by renormalization condition. In the case of $\Gamma=\gamma_{t}$, the above definition is simply

[TABLE]

which is similar to a discrete scale-dependent anomalous dimension $\partial\log\left(Z_{\gamma_{t}\gamma_{t}}(z,p)\right)/\partial p$. Through the dependence of $\zeta$ on $p^{R}$ slightly away from $p$, we can understand how well the leading order perturbation theory is able to describe the exact non-perturbative determination on the lattice. It is important to stress that apart from understanding non-perturbative renormalization of qPDF in this way, we are also essentially comparing one of the steps in the LaMET formalism that is calculable on lattice. Hence, any agreement/disagreement we observe quantifies the limitations of the leading order LaMET. In perturbation theory, $\zeta$ is the ratio of the one-loop perturbative correction to $q(z,p)$ to its tree-level value. This expression for $\zeta$ has been calculated, and it is given by 111The formula differs from the one given in Stewart and Zhao (2018) due to the issue of order of $\epsilon=0$ limit in dimensional regularization and the $z=0$ limit. We thank Yong Zhao for communicating the corrected result to us.

[TABLE]

where $H(x,p)$ is the 1-loop correction term to the bare qPDF, and the two terms in the right hand side come from the bare and RI-MOM renormalization counter terms respectively 222The function $H(x,p)$ here is referred to as $h(x,p)$ in Ref Stewart and Zhao (2018). We reserve $h$ to refer to qPDF matrix element as is the convention.. The functional forms of $H(x,p)$ for $\gamma_{z}$ and $\gamma_{t}$ isovector qPDFs are given in Stewart and Zhao (2018); Liu et al. (2018) and therefore we do not provide them here. The asymptotic $3/(2|x|)$ behavior of the bare and the RI-MOM counter term, that contributes to the UV divergence when integrated over $x$, gets exactly canceled and we obtain a UV finite and renormalized result for $\zeta$. In the discussions below, we will consider the cases with $p_{z}=p_{z}^{R}$ and $p_{z}\neq p_{z}^{R}$ separately. In the above leading order formula, the scale at which $\alpha_{s}$ has to be evaluated is arbitrary. Therefore, we vary $\alpha_{s}$ by changing the scale from $0.5p_{z}^{R}$ to $2p_{z}^{R}$ through the 1-loop running, and quote this variation as an uncertainty in the perturbative results below. On the lattice side, we determine $\zeta(z,p,p^{R})$ using the non-perturbatively determined $Z$-factors. In order to estimate the lattice spacing effects, we determined $\zeta$ using two different lattice spacings; $a=0.04$ fm is shown as filled symbols and $a=0.06$ fm is shown as open symbols in the various plots that follow.

In the left and right panels of Fig. 10, we show the typical dependence of $\zeta_{\gamma_{t}}(z)$ and $\zeta_{\gamma_{z}}(z)$ respectively, as a function of $z$ when $p_{\perp}^{R}$ differs slightly from the transverse quark momentum $p_{\perp}$, while the longitudinal components $p_{z}$ and $p_{z}^{R}$ are the same. Using Eq. (29), we calculated the prediction from leading order perturbation theory for $\zeta(z)$ at the same values of momenta. The uncertainty bands for the perturbative result are shown in Fig. 10 along with the actual lattice data at the two different lattice spacings that are shown using symbols. For the data shown in Fig. 10, the longitudinal components $p_{z}$ for the two lattice spacings are exactly 1.92 GeV, but the transverse components $p_{\perp}$ are only approximately the same between the two due to the constraints of allowed momenta on the two different lattice volumes i.e., $p_{\perp}=1.49$ GeV for $a=0.06$ fm and $p_{\perp}=1.67$ GeV for $a=0.04$ fm. To take care of this slight offset in $p_{\perp}$ between the two lattice spacings, we have distinguished the perturbative results corresponding to $a=0.04$ fm as bands enclosed by solid lines, and similarly for $a=0.06$ fm as bands enclosed by dashed lines. It can be seen that the two perturbative results are not very sensitive to this difference in $p_{\perp}$ assuring us that whatever change we observe between the data at two different $a$ is mainly due to the change in $a$. We observe from the plots that the leading order perturbation theory captures the qualitative $z$-dependence of both the real and imaginary parts of $\zeta$ when $p_{\perp}^{R}$ is changed from $p_{\perp}$. Surprisingly, the 1-loop result seems to work better for $\zeta_{\gamma_{z}}$ than for $\zeta_{\gamma_{t}}$. In the case of $\gamma_{t}$, one can certainly see a large lattice spacing effect with the movement of data towards the 1-loop result as the lattice spacing is reduced, while in the case of $\gamma_{z}$, one can already see a consistency with the one loop result at the lattice spacings that we use. Thus, it opens up a question on whether the $\gamma_{z}$ qPDF fares worse compared to the $\gamma_{t}$ qPDF simply due to the presence of mixing with the scalar or whether $\gamma_{z}$ qPDF might eventually show better perturbative convergence and lesser lattice spacing dependence in spite of its other disadvantages.

In Fig. 11, we concentrate on the renormalization flow of $\zeta_{\Gamma}$ at fixed small value of $z=0.12$ fm. The two panels show the dependence of $\zeta_{\gamma_{t}}$ and $\zeta_{\gamma_{z}}$ on $p^{R}_{\perp}$ which is changed around $p_{\perp}$. As before, we keep $p_{z}=p^{R}_{z}=1.92$ GeV. The 1-loop result is able to capture the qualitative trend of the flow in both $\gamma_{z}$ and $\gamma_{t}$. For both the cases, we can see that the reduction of lattice spacing leads to a better agreement with the 1-loop result. Having discussed the cases where $p_{z}=p_{z}^{R}$, we now study the dependence of $\zeta$ on $p_{z}^{R}\neq p_{z}$, while keeping $p_{\perp}=p^{R}_{\perp}$. We show the $z$-dependence of $\zeta_{\gamma_{z}}$ when $p_{z}^{R}=1.5p_{z}$ in Fig. 12. We find the perturbative result to have the same qualitative behavior as the lattice data. Putting together the various observations in this section, we found only an overall qualitative agreement between the lattice results on $\zeta$ and the one-loop perturbative results. When the lattice spacing is reduced, we found the agreement to get better. It remains to be seen what the effect of including higher order corrections in the perturbative result for $\zeta$.

V.2 A way to classify quark-antiquark separations as perturbative or

nonperturbative

Physically, one would expect that for a well-separated quark-antiquark with $z>1$ fm, one would start seeing traces of nonperturbative physics in the qPDF. Quantifying the advent of non-perturbative physics for large enough $z$ at finite quark/hadron momentum is important with regard to the extraction of PDF since the real-space qPDF at all $z$ enter the computation of its Fourier transform. A simple first approximation to study this effect is the following. In free theory, the qPDF with external quark states is a pure wave $e^{ip_{z}z}$. We expect, to a first approximation, that the effect of nonperturbative physics is to damp this pure wave via an inverse screening length $m_{\rm scr}\sim{\cal O}(\Lambda_{\rm QCD})$. Thus, we model the bare quark qPDF as

[TABLE]

where we have removed the UV divergent piece $e^{-c|z|}$ from the qPDF and defined the left-over exponent $m_{\rm scr}$ as a physical scale. We have also accounted for $\omega\neq p_{z}$ in the interacting theory since the quark can lose momentum by emitting gluons. There could be remnant non-trivial dependence of the amplitude $A$ on $z$, which we assume to be sub-leading compared to the leading damped oscillatory behavior and ignore them in the discussion here. There is an ambiguity in $m_{\rm scr}$ depending on the scheme used to determine the divergent piece $c$. Since the values of $c$ determined from the static quark potential method ensures that the renormalization factors after the removal of $e^{-c|z|}$ are ${\cal O}(1)$ at smaller $z$ in Fig. 8, the choice of separation of the exponential suppression factor into a divergent and physical scales as defined in Eq. (30) is well motivated in this scheme. In Fig. 13, we show the bare quark qPDF $q_{\gamma_{t}}(z,p)$ for $p_{\perp}=1.49$ GeV and $p_{z}=1.29$ GeV determined on the $a=0.06$ fm ensemble. The short distance can simply be described by a pure oscillatory $e^{i\omega z}$ behavior which is shown using the dashed curves (with $\omega=0.85p_{z}$ for the case shown). The solid curves in the figure correspond to the ansatz in Eq. (30) which describes the data at larger $|z|$ well. Without dwelling further on finding the best parametrization of the lattice data that asymptotically behaves like Eq. (30), we simply define an effective $z$-dependent $\omega$ and $m_{\rm scr}$ through

[TABLE]

In Fig. 14, we show the behavior of $m_{\rm scr}$ and $\omega$ as a function of $z$ as extracted from $q_{\gamma_{t}}(z,p)$. We have chosen different set of $p_{z}$ and $p_{\perp}$ to show the dependence on $p_{z}$ at fixed $p_{\perp}$ and vice versa. From the top panel, we see that $\omega/p_{z}$ is below 1 for $z<0.4$ fm and seems to approach a plateau closer to 1 for $z>0.4$ fm. While the values of $\omega$ at short distances depend on $p_{z}$ and $p_{\perp}$, the approach to $\omega\approx p_{z}$ is universal. We observed this behavior when we used $q_{\gamma_{z}}(z,p)$ as well. A physical reasoning for this observation could be that at shorter $z$, the quark has the ability to radiate a gluon, and at distances $z>0.4$ fm there is effectively a dressed quark carrying all the momentum. In the bottom panel of Fig. 14, we have shown the effective screening mass $m_{\rm scr}$. In the plots, we have only shown the data where 1-HYP smeared Wilson line was used. For this case, we subtracted $ca=0.1586$ in Eq. (31) to get $m_{\rm scr}$. One can clearly see the emergence of non-zero $m_{\rm scr}\approx 300$ MeV for $|z|>0.5$ fm which is in the typical ${\Lambda}_{\rm QCD}$ scale. When we repeated this using quark qPDF with unsmeared Wilson line, we found the results to be consistent with the data shown in Fig. 14 after we subtracted out $ca=0.3687$ corresponding to unsmeared Wilson line. This assures us that the observed $m_{\rm scr}\approx 300$ MeV is a real physical scale independent of the self-energy divergence of the Wilson line. This signals the significant presence of a confinement scale beyond $z\approx 0.5$ fm. Also, the near plateauing of both $\omega$ and $m_{\rm scr}$ for these larger $z$ indicates that a simple physically motivated ansatz in Eq. (31) offers a surprisingly good description of the actual non-perturbative data. One could have expected this simply from observing the large $|z|$ part of Fig. 13. It remains to be seen if this observation can be used advantageously in improving the LaMET matching at finite moderately large $p_{z}$.

VI From renormalized quasi-PDF to PDF

VI.1 On obtaining the valence PDF using isovector $u-d$ qPDF of pion

Having determined the renormalized qPDF we can now discuss the matching between qPDF and PDF as well as the determination of pion PDF from our lattice results. We computed the $u-d$ qPDF matrix element of a pion which in practice we obtained from the real part of the connected piece of the $u$ quark qPDF matrix element. Now, we discuss how the $u-d$ qPDF and PDF are related to the valence PDF of pion.

The $u$ and $d$ quark distributions, $f_{u}(x)$ and $f_{d}(x)$, as determined using Eq. (1) has support from $x=-1$ to 1. One can make connection with the conventional, separately defined quark distributions $Q_{u,d}(x)$ and the anti-quark distributions $Q_{\bar{u},\bar{d}}(x)$ that are non-zero only between $x=0$ and 1, through the relation

[TABLE]

Therefore, $f_{u,d}(x)$ contains information on both the quark as well as the antiquark distributions in the positive and negative regions of $x$ respectively. Let us first focus on $x\geq 0$. In the isospin symmetric case we are considering, $Q_{u}(x)=Q_{\bar{d}}(x)$ and $Q_{\bar{u}}(x)=Q_{d}(x)$. Therefore, for the positively charged pion $Q_{u}(x)-Q_{d}(x)=Q_{u}(x)-Q_{\bar{u}}(x)=f^{\pi,u}_{v}(x)$ is the valence u-quark distribution. Again, due to the isospin symmetry, the $u$ and $\bar{d}$ valence distributions are the same as $f^{\pi,u}_{v}(x)=f^{\pi,{\bar{d}}}_{v}(x)=f^{\pi}_{v}(x)$. However, unlike the valence quark distribution, the isotriplet $u-d$ PDF $f_{u-d}=f_{u}(x)-f_{d}(x)$ satisfies $f_{u-d}(|x|)=f_{u-d}(-|x|)$ and it has support from $-1$ to 1. That is,

[TABLE]

Therefore, one can obtain the $u-d$ quark distribution, and from it, one can obtain $f_{v}^{\pi}(x)$ from $x\in[0,1]$, or equivalently, from $[-1,0]$.

By applying the matching formula on $f_{u}(x)$ and $f_{d}(x)$ separately and taking the difference to obtain the $u-d$ RI-MOM qPDF, we now try to learn what is expected for this qPDF. Writing down only the $x/y$ dependence for the sake of brevity and keeping the dependence on $yP_{z}/P_{z}^{R}$, $(P^{R}/P^{R}_{z})^{2}$ and factorization scale $\mu$ implicit, the one-loop contribution to the matching kernel $C(x/y)$ from the RI-MOM to ${\overline{\rm MS}}$ scheme consists of two terms: ${\mathcal{F}}_{1}(x/y)$ and ${\mathcal{F}}_{2}(1+\eta^{\prime}(x-y))$ with $\eta^{\prime}=P_{z}/P_{z}^{R}$. The expressions 333In Stewart and Zhao (2018), the terms ${\mathcal{F}}_{1}$ and ${\mathcal{F}}_{2}$ are referred to as $f_{1}$ and $f_{2}$, respectively for ${\mathcal{F}}_{1,2}$ depend on the choice of $\Gamma$ ($\gamma_{z}$ or $\gamma_{t}$) Liu et al. (2018); Stewart and Zhao (2018). Furthermore, ${\mathcal{F}}_{2}$ depends on the projection method of the RI-MOM scheme. Using the matching formula Stewart and Zhao (2018); Liu et al. (2018) on $f_{u}(x)$ and $f_{d}(x)$ to obtain the qPDFs $\tilde{q}_{u}(x)$ and $\tilde{q}_{d}(x)$,

[TABLE]

The above equation includes both the sea and valence quarks, and there will be mixing with the gluon PDF which are included in the “$\cdots$” part. In the above convolution, the vector current conservation is ensured by the plus function defined as

[TABLE]

such that any extra variable that ${\mathcal{F}}_{1,2}$ will depend on are held fixed in the above integral. Since the matching between qPDF and PDF is linear, the $\tilde{q}_{u-d}=\tilde{q}_{u}(x)-\tilde{q}_{d}(x)$ is simply obtained as

[TABLE]

with the terms in “$\cdots$” in Eq. (36) exactly canceled between the $u$ and $d$ terms. This is the matching relation we use to obtain the $u-d$ PDF from $u-d$ qPDF. Using the $u-d$ PDF, we obtained the valence PDF as discussed above. While $f_{u-d}(|x|)=f_{u-d}(-|x|)$, it is also true that $\tilde{q}_{u-d}(|x|)\neq\tilde{q}_{u-d}(-|x|)$ in the RI-MOM scheme. One way to understand this is from the fact that the bare qPDF matrix element is purely real while the RI-MOM renormalization factor is in general complex, thereby making the renormalized qPDF matrix element complex. One can see this by starting from the matching convolution above and find that,

[TABLE]

is non-zero due to an RI-MOM specific term ${\mathcal{F}}_{2}$, while the terms containing ${\mathcal{F}}_{1}$ cancel due to their dependence only on $|P_{z}|$. In other schemes such as the ${\overline{\rm MS}}$, this symmetry about $x=0$ would be preserved by matching because the corresponding factorization formulas depend on renormalization/regularization scales through combinations such as $\mu^{2}z^{2}$ Izubuchi et al. (2018). In the RI-MOM scheme there are two renormalization scales, $P^{R}$ and $P_{z}^{R}$, and since the $z$-direction is special the above statement does not hold. Thus, it is important to capture this asymmetry in the qPDF, or equivalently to describe both the real and imaginary parts of the RI-MOM renormalized pion qPDF from matching. We use the matching kernel corresponding to the $\not{p}$-projection in the results to be discussed next.

VI.2 Numerical results on pion valence PDF from matching

The one-loop perturbative matching relates the Fourier transform $\tilde{q}(x,P_{z},P^{R})$ of the renormalized RI-MOM real-space qPDF matrix element $h^{R}(z,P_{z},P^{R})$, and the ${\overline{\rm MS}}$ PDF $f(x,\mu)$ at factorization scale $\mu$. The relation is through the convolution in Eq. (3). There are two approaches to consider here:

One can parametrize the real space data $h^{R}(z,P_{z},P^{R})$ over the range $z$ where one has the lattice data and then model the dependence of $h^{R}(z,P_{z},P^{R})$ over $z$ extending to infinity where data does not exist (c.f., Chen et al. (2017); Lin et al. (2018b)). Using such a parametrization, one can obtain its Fourier transform $\tilde{q}(x,P_{z},P^{R})$. Since, the matching is only up to ${\cal O}({\alpha_{s}})$, one can invert the relation Eq. (3) by replacing $f\leftrightarrow\tilde{q}$ and $\alpha_{s}\to-\alpha_{s}$. Thereby, one can obtain $f(x,\mu)$. In this method, one does not control what values of $z$ enter the Fourier transform and one could question the validity of perturbation theory for $z>1$ fm. 2. 2.

One can start from a phenomenologically motivated $n$-parameter family of PDFs $f(x,\mu;a_{1},\ldots a_{n})$. Through Eq. (3), one can obtain qPDF $\tilde{q}(x,\mu;a_{1},\ldots a_{n})$, and thereby, obtain a family of real space space qPDF matrix elements $h^{R}(z,P_{z},P^{R};a_{1}\ldots a_{n})$. Using this, one can fit the parameters $(a_{1},\ldots,a_{n})$ so as to best describe the real space lattice data over a range $z$. This method was used in the case of lattice cross-section approach in Sufian et al. (2019). Since the model PDFs are not predictions from QCD, the model dependence enters the analysis and one has to rely on the prior that experimentally determined PDFs are indeed very well described by such a family of PDFs. However, the advantage of this method is that one can precisely control the range of $z$ that enters the analysis, and one also does not have to invert the matching convolution.

From our observation on how the 1-loop perturbation theory fails to capture the quark qPDF quantitatively even at short distances and from the observation of significant nonperturbative screening effects beyond $z=1$ fm, we think it is important to be in control of what values of $z$ enter the convolution and hence, in this paper we take the second approach. Also, due to the loss of signal to noise ratio for $z>1$ fm, we found Fourier transforming the noisy data to be challenging without introducing unwanted wiggles in $\tilde{q}(x)$ at larger $x$.

To be on par with the experimental extraction of PDFs, one should use sophisticated methods such as the usage of neural networks to choose the set of model PDFs to start with (c.f., Karpie et al. (2019)). We defer such an analysis to a future work and instead, we use a simple two-parameter phenomenologically motivated functional form for the valence PDF:

[TABLE]

for $x\in[0,1]$ and zero elsewhere. As we will see below, such a form is enough to describe our lattice data. One can fix the coefficient $A$ through a stringent condition $\int_{0}^{1}f^{\pi}_{v}(x)dx=1$. Instead, we use a more conservative constraint on $A$ using $\int_{0}^{1}f^{\pi}_{v}(x)dx=h^{R}(z=0,P_{z},P^{R})$ to allow for sample by sample fluctuations in $h^{R}(z=0,P_{z},P^{R})$ close to 1 and fold this into the error estimate. It should be noted that the valence PDF of pion determined from the experimental data by the JAM collaboration Barry et al. (2018) can be well described by such a two parameter ansatz, for example with $a=-0.407$ and $b=1.12$ at $\mu=3.2$ GeV.

Using the above valence PDF, we construct the $u-d$ PDF as

[TABLE]

with $x\in[-1,1]$ and zero elsewhere. Through the convolution of $f_{u-d}(x)$ with the matching kernel, we obtain $\tilde{q}_{u-d}(x;a,b)$, which in turn we use to construct the real space qPDFs $h(z;a,b)=\int_{-\infty}^{\infty}\tilde{q}(x;a,b)e^{ixP_{z}z}dx$. We will refer to these functions $h(z;a,b)$ as the two-parameter family of phenomenologically motivated qPDF matrix elements. With the set of $h(z;a,b)$ from a range of $a$ and $b$, we can fit the parameters $a$ and $b$ to the data by minimizing either $\chi^{2}_{r}$ or $\chi^{2}_{ri}$ below:

[TABLE]

In the above equations, $[-z_{\rm max},z_{\rm max}]$ specifies the fit range. The statistical errors on the real and imaginary parts of the lattice data $h^{R}(z)$ is $\sigma_{r}(z)$ and $\sigma_{i}(z)$ respectively. To account for any systematic errors coming from the higher order corrections in $\alpha_{s}$ in the matching kernel, we determine $h(z;a,b)$ from $f^{\pi}_{v}(x;a,b)$ by varying the value of $\alpha_{s}$ in the matching kernel from $\alpha_{s}(\mu/2)$ to $\alpha_{s}(2\mu)$ though the 1-loop running. The corresponding changes in the real and imaginary parts of $h(z;a,b)$ are denoted as $\sigma_{r}^{\rm pert}(z)$ and $\sigma_{i}^{\rm pert}(z)$ respectively, and we include these uncertainties in the matched result in the $\chi^{2}$. If the matching was exact, then by fitting only the real part by minimizing $\chi^{2}_{r}$ would automatically guarantee that the imaginary part also agree with the data. Therefore at any finite order matching, the fits obtained by minimizing $\chi^{2}_{r}$ and $\chi^{2}_{ri}$ will in general be different. For the results shown below, we used $\chi_{ri}^{2}$ in order to obtain the PDF that best describes both the real and imaginary parts of the real space qPDF, but we also used $\chi^{2}_{r}$ and found it to lead to consistent results, but with larger uncertainties. We did not include the correlations between the data at different $z$ for the primary reason that it is difficult to keep these correlations intact in the process of excited state extrapolations. It also helps us to easily incorporate the effect of $\sigma^{\rm pert}$ from non-statistical origin in the analysis, and in treating the real and imaginary parts of the renormalized matrix elements as two distinct pieces of data as is the case in the context of matching. We determined the errors on the fit parameters through the bootstrap analysis.

In Fig. 15, we show the fitting procedure for $\gamma_{t}$ qPDF. In the top panels, we show the $P_{z}=0.86$, 1.29 and 1.72 GeV real-space RI-MOM pion qPDF matrix elements from left to right. The symbols are the actual lattice data. The solid and patterned red (blue) bands are 1-$\sigma$ error-bands of the real (imaginary) parts of the fitted real space qPDF matrix element that best fits the data over the range $[-z_{\rm max},z_{\rm max}]$ for $z_{\rm max}=1.44$ fm and 0.72 fm respectively. The agreement with both the real and imaginary parts of the lattice data is noteworthy. In fact, we find the qPDF matrix element as inferred from the JAM PDF Barry et al. (2018) is able to explain the lattice data well for the entire range of $z$ at the two largest momentum. In the bottom panels, we show the process leading from model PDF to the real space qPDF matrix elements shown in the top panels. In order to avoid cluttering the figure, we have shown only the mean value of $f_{u-d}(x)$ (shown as dashed lines) while we have shown the error bands for the qPDF $\tilde{q}_{u-d}(x)$ as obtained through the 1-loop matching. The colors red and blue in the bottom panels correspond to the fits with $z_{\rm max}=0.72$ and 1.44 in the top panels respectively. As one can see, we started from a symmetric $u-d$ PDF by construction and matching introduces an $x\to-x$ asymmetry. After Fourier transformation, this asymmetry leads to the imaginary part in the real space data in the top panels which captures the lattice data to a good accuracy. For both the real-space as well as in $x$ space, we find no significant difference between using $z_{\rm max}=1.44$ fm and $0.72$ fm in the fits. We could infer that within the precision of our numerical results, the non-perturbative effects at $z\approx 1$ fm that we found using quark qPDFs is not important. Therefore, we show results for an intermediate $z_{\rm max}=0.98$ fm in the results below. When we repeated this analysis by minimizing $\chi^{2}_{r}$, we found the estimates to be consistent with the above, but with larger uncertainties.

In Fig. 16, we show our results for $f_{v}^{\pi}(x,\mu)$ and $xf_{v}^{\pi}(x,\mu)$ at the factorization scale $\mu=3.2$ GeV using the procedure described above at our two largest pion momenta $P_{z}=1.29$ GeV and $1.72$ GeV starting from $\gamma_{t}$ qPDF. For each case, we overlay the results from two different RI-MOM scales $P^{R}$ in order to show the scatter as a systematic error in our estimates. We find the $P^{R}$ dependence to be minor compared to the error bands (we repeated the analysis with multiple other values of $P^{R}$ that are not shown and only minor scatter with respect to $P^{R}$ was seen). We also show the result from the JAM collaboration Barry et al. (2018) for the pion valence PDF at the same factorization scale as the black solid line, which lies within the statistical and systematic uncertainties of our estimates. In the left panels showing $f_{v}^{\pi}(x)$, this overall agreement can be seen even up to smaller $x$, but one has to be cautious of this agreement for $x\lessapprox\Lambda_{\rm QCD}/P_{z}\approx 0.2$ for the two highest pion momenta we use. By construction, in our fitting procedure $f_{v}^{\pi}(x)$ has support only from 0 to 1 without any necessity to recover this condition in the infinite $P_{z}$ limit. However, the values of exponent $b$ closer to zero are also allowed thereby leading to a wider error band closer to $x$=1. This seems to be consistent with the observation in Ref Chen et al. (2018b) that the PDF obtained from qPDF through the inverse one-loop matching (approach-1) vanishes at about $x\approx 1.2$. We see our $P_{z}=1.29$ GeV and 1.72 GeV estimates to be consistent albeit with a significant increase in error at the largest momentum.

In Fig. 17, we summarize the information in Fig. 16 by showing the 1-$\sigma$ ellipses (whose $x$ and $y$ projections give the marginal 68% confidence intervals of the exponents $a$ and $b$ respectively). In this figure, the dashed and continuous ellipses are for $P_{z}=1.29$ and 1.72 GeV respectively. The ellipses for different $P^{R}$ are distinguished by the colors, with the color code being the same as in Fig. 16. The $P_{z}=1.29$ GeV data offers a stronger constraint on allowed region of $(a,b)$ than the noisier $P_{z}=1.72$ GeV. In this plot, the JAM estimate is the black point. The JAM data is well within the $P_{z}=1.72$ GeV ellipses while the $P_{z}=1.29$ GeV data seems to favor slightly smaller exponent $b$. However, these differences are well within $2\sigma$. Even though our lattice data has large errors on the exponents $a$ and $b$ individually, the data offers a tight constraint on the combined allowed region. In particular, the principal component of this correlation between $a$ and $b$ points directly at the JAM data implying that if one fixes the exponent $a$ to be from the experiment, then the best value of $b$ would also be closer to that from the experiment. To understand this better, we have also shown the line of constant value of first moment of the valence PDF, $\langle x\rangle=\int_{0}^{1}xq_{v}^{\pi}(x)dx$, set to 0.215 as inferred from the JAM data. It is clear that the 1$\sigma$ ellipses are oriented along this line, which means that qPDF determines $\langle x\rangle$ robustly and this in turn provides a strong constraint in the allowed PDFs. Not surprisingly, we do find consistent values of $\langle x\rangle=$ 0.21(2) and 0.22(3) from the $P_{z}=1.29$ and 1.72 GeV estimates. It should be noted that the moments of pion PDF have also been directly determined without the usage of LaMET formalism Best et al. (1997); Guagnelli et al. (2005); Capitani et al. (2006); Abdel-Rehim et al. (2015); Oehm et al. (2019) and similar values for the first moment for the pion were obtained, but at slightly different values of $\mu^{2}$ than used here.

The exponents $a$ and $b$ were also recently obtained using the lattice cross-section approach Sufian et al. (2019) which used current-current correlators, with the matching implemented at tree-level. Here, the exponents were estimated as $a=-0.34(31)$ and $b=1.93(68)$ which are consistent with the region allowed at the largest momentum in Fig. 17. It is worth noting that there are indications from next-to-leading-logarithmic soft gluon resummation calculation Aicher et al. (2010), the Dyson-Schwinger equation Nguyen et al. (2011); Chen et al. (2016); Bednar et al. (2018) and light-front holographic QCD de Teramond et al. (2018) that the value of exponent $b$ could be approximately 2 as expected from perturbative counting rule (c.f., Brodsky (2005)), whereas a chiral quark model analysis Ruiz Arriola (2002); Broniowski and Ruiz Arriola (2017) suggests a value of $b$ closer to 1. It will be interesting to see if a similar implementation of an improved matching kernel could lead to a softer large $x$ behavior for the pion than what is observed using the 1-loop qPDF matching here and perhaps in Chen et al. (2018b). In fact, a general consideration of power corrections to qPDF Braun et al. (2019) revealed the presence of the form $\Lambda^{2}_{\rm QCD}/\left((1-x)x^{2}P^{2}_{z}\right)$ implying higher values of $P_{z}$ might be required in order to correctly describe physics close to $x=1$, and this might be the effect which we are finding. Similar conclusions have also been obtained in 2d QCD Ji et al. (2019).

Due to the larger errors in the $\gamma_{z}$ qPDF attributed mostly to the steep excited state extrapolations, we use the $\gamma_{z}$ qPDF to provide a consistency check of our calculations instead. For this, we use our above best estimates of the PDF obtained using the $\gamma_{t}$ qPDF to get the corresponding prediction for the real space $\gamma_{z}$ qPDF matrix element through a convolution with the appropriate matching kernel. In the top and bottom panels of Fig. 18, we show such a comparison between the actual real space data of $\gamma_{z}$ qPDF (data points) along with the prediction from our estimated PDF (bands) for pion momenta $P_{z}=1.29$ and 1.72 GeV. We find good descriptions of the real part of the RI-MOM $\gamma_{z}$ qPDF at both the pion momenta with a slight tension between the imaginary parts. From our discussion on the excited state contamination, it is important to first gain better control of the larger excited state contamination in $\gamma_{z}$ qPDF before one can investigate the effect of one-loop matching on this rather small discrepancy.

VII Conclusions

We studied pion PDF in the framework of LaMET, which relates the qPDF to PDF through the matching convolution in Eq. (3). For this, we used a small lattice spacing $a=0.06$ fm. We carefully examined the effects of excited states using the two- and three-state exponential fits of the relevant 2-point and 3-point functions as well as the summation method. For our final analysis we used two momenta $P_{z}=1.29$ GeV and $P_{z}=1.72$ GeV. We found the qPDF defined using $\Gamma=\gamma_{t}$ was better determined compared to the $\gamma_{z}$ qPDF in the lattice calculation for the following reasons: smaller statistical error, relatively smaller excited state extrapolation leading to more robust result for the matrix element as well as due to the absence of mixing. Therefore, we focused on the analysis of the $\gamma_{t}$ matrix element.

The pion qPDF was non-perturbatively renormalized using the RI-MOM scheme by calculating the matrix elements of qPDF operator with off-shell quark states in the Landau gauge for different separations $z$. For these calculations we also used finer lattices with the lattice spacing $a=0.04$fm. We performed the comparison of this matrix element in the Landau gauge with 1-loop perturbative calculations in the RI-MOM scheme and found a qualitative agreement for $z<0.3$ fm. For the smaller lattice spacings, $a=0.04$ fm we even found quantitative agreement with the 1-loop result for sufficiently small $z$. We also explored the role of non-perturbative effects in the calculation of the off-shell matrix element. The real part of the RI-MOM renormalization coefficient is close to one, while the imaginary part is close to zero once the divergent self energy part of the Wilson line is removed. We pointed out that the RI-MOM renormalization procedure leads to an asymmetry in the iso-vector pion qPDF $\tilde{q}(x,P_{z},p_{z}^{R},\mu_{R})$ around $x=0$, while other renormalization procedures lead to qPDF that is symmetric around $x=0$.

From the renormalized qPDF, we determined the valence quark pion PDF using the 1-loop perturbative matching of the $\gamma_{t}$ qPDF, which we implemented through a fit to the phenomenologically motivated $x^{a}(1-x)^{b}$ functional form for the valence PDF. We found our results for the pion valence PDF using the two largest pion momenta were consistent with each other, though the statistical errors are rather large. An overall agreement with the results obtained recently by the JAM collaboration Barry et al. (2018) was seen. We found our result for the PDF to capture the first moment $\langle x\rangle$ more robustly than the small-$x$ and large-$x$ exponents, $a$ and $b$ themselves. We used the $\gamma_{z}$ qPDF matrix elements to provide an internal consistency check by comparing to the expectation from our estimates of the PDF and a satisfactory agreement was seen. From our analysis it is clear that the dominant source of errors in the PDF determination is the statistical error of the lattice calculations. It will be necessary to significantly increase the statistics in the future lattice calculations. Future high statistics lattice calculations will be important for an accurate determination of the pion PDF as well as testing of the LaMET approach around small $x$.

Acknowledgments

We thank P.C. Barry for sharing the JAM valence PDF data, and J. H. Weber for providing us the Wilson line data. We thank R. S. Sufian and Y. Zhao for fruitful discussions. This material is based upon work supported by: (i) The U.S. Department of Energy, Office of Science, Office of Nuclear Physics and High Energy Physics through the Contract No. DE-SC0012704; (ii) The U.S. Department of Energy, Office of Science, Office of Nuclear Physics and Office of Advanced Scientific Computing Research within the framework of Scientific Discovery through Advance Computing (ScIDAC) award Computing the Properties of Matter with Leadership Computing Resources; (iii) The Brookhaven National Laboratory’s Laboratory Directed Research and Development (LDRD) project No. 16-37. S.S. is supported by the National Science Foundation under CAREER Award PHY-1847893. This research used awards of computer time provided by: (i) The INCITE program at Oak Ridge Leadership Computing Facility, a DOE Office of Science User Facility operated under Contract No. DE-AC05- 00OR22725; (ii) The USQCD consortium at its Brookhaven National Laboratory and Jefferson Laboratory computing facilities. The gauge configurations used in this study were generated using awards of computing time provided by: (i) The INCITE program at Argonne Leadership Computing Facility, a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-06CH11357.

Appendix A Coulomb-gauge Gaussian and momentum (boosted) smearing

In order to create hadron interpolating operators that have a good overlap with the corresponding ground states, quark field smearing is typically required. The amount of applied smearing is tuned to produce spatial quark distributions of roughly of the same spatial size as the hadron. Gauge-covariant Wuppertal (Gaussian) smearing Gusken (1990) is commonly used for this purpose. However, calculation of quasi- and pseudo-PDFs and high-momentum hadron structure and spectrum in general requires lattices with small lattice spacing. Keeping the physical size of smeared quark distributions the same becomes a numerical challenge on finer lattices because it requires larger numbers of smearing iterations. For this reason, we use Gaussian shape smearing in a fixed (Coulomb) gauge that can be performed efficiently through a convolution with a Gaussian profile kernel,

[TABLE]

In the free-field case, this kernel corresponds to the Wuppertal smearing operator $(1+\frac{w^{2}}{4N}\Delta_{\text{sp}})^{N}$, where $\Delta_{\text{sp}}$ is the spatial Laplacian and $w^{2}=2w_{\text{CG}}^{2}$. The value for the width $w_{\text{CG}}$ is chosen to match the mean-squared radius $\langle r^{2}\rangle=3w_{\text{CG}}^{2}$ to that of the optimal Wuppertal-smeared quark sources. First, we fix the Coulomb gauge

[TABLE]

where $\Omega^{C}_{x}$ is the gauge transformation to the Coulomb gauge, which minimizes the functional

[TABLE]

(for the Coulomb gauge, $\mu\neq t$ and the functional is minimized independently on each time slice). The numerical implementation is identical to the algorithm used for fixing the Landau gauge in NPR calculations. Application of the smearing kernel requires two 3D Fourier transformations

[TABLE]

which is accelerated with offloading matrix-matrix products to GPU.

Incorporating momentum (boosted) into Coulomb-gauge Gaussian smearing amounts to translation of the kernel in the momentum space,

[TABLE]

In a periodic finite volume, care must be taken to avoid spatial discontinuities in the boosted smearing kernel (Eq. (52)). Such discontinuities may arise because the optimal boosted smearing momentum $\vec{k}$ typically does not conform to finite-volume momentum quantization $\vec{k}=2\pi\frac{\vec{n}}{\vec{L}}$ and the phase factors $e^{i\vec{k}\vec{x}},\,e^{i\vec{k}\vec{y}}$ do not satisfy periodic boundary conditions. The solution is to define the smearing kernel in the momentum space as ${\mathcal{S}}^{(\vec{k})}_{\vec{x},\vec{y}}=\sum_{\vec{p}}e^{i\vec{p}(\vec{x}-\vec{y})}\,e^{-\frac{1}{2}w_{CG}^{2}(\vec{p}-\vec{k})^{2}}$, where the momentum difference $(\vec{p}-\vec{k})$ is understood as the shortest distance between $\vec{p}$ and $\vec{k}$ in the Brillouin zone. Such choice leads to a smooth distribution in the momentum space and respectively smooth and continuous smearing kernel in the coordinate space.

Finally, it is important to note that the smearing kernel in Eq. (52) is Hermitian (as an operator acting in the [coordinate $\otimes$ color] space),

[TABLE]

which is similar to the (boosted) Wuppertal smearing operator and important for computing symmetric hadron correlation functions.

Appendix B Meson correlation functions with boosting

We use the interpolating operator for the $\pi^{+}=\bar{d}u$ meson

[TABLE]

which is constructed from smeared quark fields

[TABLE]

where the spinor matrix $\Gamma_{M}=\gamma_{5}$. The Hermitian-conjugated (creation) meson operator is

[TABLE]

where $\overline{\Gamma}_{M}=\gamma_{4}\Gamma_{M}\gamma_{4}=(-\gamma_{5})$. The meson two-point correlation function with boost-smeared source and sink and momentum projection at the sink is444 Unless explicit summation is performed, implicit summation over repeated coordinate indices $x^{\prime},x^{\prime\prime},y^{\prime},y^{\prime\prime}$ is assumed, as well as all over omitted spin and color indices.

[TABLE]

where $Q^{q}_{x,y}=\langle q_{x}\bar{q}_{y}\rangle$ and $\tilde{Q}^{q,(\pm\vec{k})}={\mathcal{S}}^{(\pm\vec{k})}Q^{q}{\mathcal{S}}^{(\pm\vec{k})}$ are unsmeared and smeared quark propagators, respectively. Note that the meson two-point function is constructed from the $u$-quark propagator $y\leftarrow x$ and the $d$-quark propagator $x\leftarrow y$ smeared with momenta $(\vec{k})$ and $(-\vec{k})$, respectively. Therefore, separate propagators for $u$ and $d$ quarks are required to construct meson correlation functions “boosted” with the total momentum $(2\vec{k})$,

[TABLE]

where “$\propto$” sign stands for additional coordinate dependence due to the boosting. The $d$-quark $x\leftarrow y$ propagator, as usual, is computed using $\gamma_{5}$-Hermiticity of the Dirac operator,

[TABLE]

where the sign of the boosting momentum is preserved due to the Hermiticity of the (boosted) smearing operator ${\mathcal{S}}$ in Eq. (53).

Repeating similar steps to for the meson three-point function with the insertion of the operator \big{[}\bar{u}\,W\,\Gamma\,u\big{]}_{z} with arbitrary $\Gamma$-matrix and Wilson line $W_{z,z+\hat{\mathcal{L}}}=(\prod_{\mathcal{L}}U)_{z,z+\hat{\mathcal{L}}}$ along path ${\mathcal{L}}$, we get

[TABLE]

where the forward propagator $F^{u}=Q^{u}\,S^{(\vec{k})}$ and the meson sink-sequential (backward) propagator $B^{\bar{d}\Gamma_{M}u(y_{4},\vec{p}^{\prime})}$ is defined as

[TABLE]

which is also computed using the $\gamma_{5}$-conjugation.

Appendix C Explicit calculation to show that

bare pion $u-d$ three point function is purely real or imaginary

In the previous appendix, we constructed the connected piece of the three point function of $\overline{u}W\Gamma u$ operator in $\pi^{+}$. If one repeats the computation using $\overline{d}W\Gamma d$ operator, one will find the disconnected piece to be the same as the one in the full $\overline{u}W\Gamma u$ three point function and hence such quark line disconnected terms will cancel in the $\overline{u}W\Gamma u-\overline{d}W\Gamma d$ isospin nonsinglet operator that we are interested in. Below, we further explain as to why only the real part of the connected $\overline{u}W\Gamma u$ three point function for $\Gamma=\gamma_{t},\gamma_{z}$ and imaginary part for $\Gamma=1$ contributes to the total isospin nonsinglet three point function. For the sake of simplicity let us take the case of point source and point sink, and take $\Gamma=\gamma_{t}$. The full expression for the $u-d$ qPDF three point function is

[TABLE]

where we do not make distinctions between $u$ and $d$ quark propagators due to isospin symmetry. Let us call the trace in first term on the right hand side as $T_{1}$ and the second trace before being conjugated as $T_{2}$. One can go from $T_{2}$ to $T_{1}$ by parity transformation $x=(\vec{x},x_{4})\to x_{p}=(-\vec{x},x_{4})$ , followed by a spatial translation $x\to x+{\mathcal{L}}$ by making use of the transformation of the Dirac propagator to be $Q_{x,y}\to\gamma_{t}Q_{x_{p},y_{p}}\gamma_{t}$ and $W_{x,x+{\mathcal{L}}}\to W^{\dagger}_{x_{p}-{\mathcal{L}},x_{p}}$ under parity. In this case, the $\gamma_{t}$ from parity transformation for $Q$ commutes with $\Gamma=\gamma_{t}$. In other cases, one should take care of the $\pm$ factor. Thus $C_{\rm 3pt}$ becomes

[TABLE]

and therefore proportional to the connected piece of $\overline{u}\Gamma Wu$, which is the first term in the above equation. We normalize the three point function such that the $u-d$ isospin charge of the pion is 1. By going through the similar calculation, one can show that the three point function is real also for $\Gamma=\gamma_{z}$ while it is is purely imaginary for $\Gamma=1$ $u-d$ pion qPDF.

Appendix D Relation between $P_{z}$ and $-P_{z}$ matrix elements

In this appendix, we derive the relation between conjugates of the matrix elements $\matrixelement{E_{n^{\prime}},P_{z}}{\mathcal{O}_{\Gamma}(z;\tau)}{E_{n},P_{z}}$ that enter the excited state contributions to the qPDF three-point function. For this, let us consider the conjugate of the simplest component of the qPDF matrix element:

[TABLE]

for $\Gamma=\gamma_{t},\gamma_{z},1$. Using the parity operator $\Pi$, the right-hand-side of the above equation becomes

[TABLE]

Defining, $\gamma_{t}\Gamma\gamma_{t}=\Phi_{\Gamma}\Gamma$ with $\Phi_{\Gamma}=\pm 1$, we have the relation

[TABLE]

with $\Phi_{\Gamma}=1$ for $\Gamma=\gamma_{t},1$ and $\Phi_{\Gamma}=-1$ for $\Gamma=\gamma_{z}$ with all $\Gamma$ being in the Minkowskian convention. Thus, we can average over $P_{z}$ and $-P_{z}$ data after taking care of appropriate factor of $\Phi_{\Gamma}$.

Appendix E Pion two point functions and energy levels

In this appendix we discuss some details of the calculations of the pion two point function. We tested several different sources for the pion. In these tests we used 50 gauge configurations. We used Gaussian sources with several steps of Wuppertal smearings as well as in Coulomb gauge (see main text). In Fig. 19 we show the effective mass for 40 and 90 steps of Wuppertal smearings as well as the Coulomb gauge Gaussian sources of size 0.3 fm. We see that 90 steps of Wuppertal smearings and Coulomb gauge Gaussian sources give similar effective masses, while the excited state contamination is larger for 40 steps of Wuppertal smearings. We also studied the two point functions for different boosted Gaussian sources, with momentum boost $k_{z}$. The corresponding effective masses are shown in Fig. 20 for $P_{z}=0.86,~{}1.29$ and $1.72$ GeV. for different values of $\zeta=k_{z}/P_{z}$. We clearly see that non-zero value of $\zeta$ improves the signal for all $P_{z}$. We also see that $\zeta=0.5$ is too small, while $\zeta=1.0$ is too large for $P_{z}=1.29$, but works well for $P_{z}=0.86$ GeV.

Appendix F Implementation of matching convolution

Here, we describe the implementation of the plus function in the matching formula such as to ensure current conservation. The matching kernel is of the form

[TABLE]

where the dependence on $P^{R}$ and $\mu$ are implicit. The first perturbative correction is a plus function that ensures the vector current conservation. The property we know of the plus-function is that $\int_{-\infty}^{\infty}dxC^{(1)}_{+}\left(\frac{x}{y},yP_{z}\right)=0$, since the second dependence of the function is independent of $x$. In order to implement the plus-function correctly, we can use the following procedure:

[TABLE]

The $x$-independent but momentum dependent coefficient $N(yP_{z})$ is

[TABLE]

where $\int_{reg}d\xi$ involves an integration over the intervals $[-\Lambda,-\epsilon]\cup[\epsilon,1-\epsilon]\cup[1+\epsilon,\Lambda]$ for some upper cut-off $\Lambda$ and a small exclusion parameter $\epsilon$. The above definition gives the usual result that

[TABLE]

with $y$ held fixed as $\Lambda\to\infty$ and $\epsilon\to 0$. The following is then true for any function $f$:

[TABLE]

leading to the vector current conservation or equivalently to the total area preservation between the qPDF and PDF. With this prescription, the matching formula becomes

[TABLE]

It is convenient to write the above formula in an explicitly vector current conservation preserving form as

[TABLE]

However, care has to be taken in the numerical regularization of the above equation to be consistent with the one in Eq. (76). That is, in the above equation, $\int_{reg}dy$ in the first integral in the right-hand side involves the range $y\in[-x/\epsilon,-x/\Lambda]\cup[x/\Lambda,x/(1+\epsilon)]\cup[x/(1-\epsilon),x/\epsilon]$ when $x>0$ and the range reversed when $x<0$. A consistent prescription for $\int_{reg}dy^{\prime}$ in the second integral in the right-hand side involves $y^{\prime}\in[-\Lambda x,-\epsilon x]\cup[\epsilon x,x(1-\epsilon]\cup[x(1+\epsilon),x\Lambda]$.

Appendix G Results on two-state extrapolations to obtain the matrix elements at all $P_{z}$

In Fig. 4 in the main text, we showed some sample results for the $t-\tau/2$ behavior and the $t\to\infty$ extrapolations of the three-point function to two-point function ratio $R(t,\tau;z,P_{z},\Gamma)$ for $\Gamma=\gamma_{t}$ and $\gamma_{z}$ at a specific intermediate value of $P_{z}=1.29$ GeV. In Fig. 21 and Fig. 22 of this appendix, we show similar results at all $P_{z}$ for $\Gamma=\gamma_{z}$ and $\gamma_{t}$ respectively using Fit(2,2).

For the case of $P_{z}=0$, special care needs to be taken. For a finite temporal extent $L_{t}$ of the lattice, ignoring the effect of periodicity due to the presence of the terms $e^{-E_{\pi}(L_{t}-t)}$ in the denominator of Eq. (10) is justified when $e^{-E_{\pi}t}\ll e^{-E_{\pi}(L_{t}-t)}$. But, one should include the effect of boundary condition if the two terms become comparable. For the largest source-sink separation $t/a=12$ we use, the contribution from the wrapping-around term, $e^{-E_{\pi}(L_{t}-t)}$, relative to $e^{-E_{\pi}t}$ for $P_{z}=0$ is 2.7%, whereas for higher $P_{z}$ it is negligible e.g., for the smallest non-zero momentum $P_{z}=0.43$ GeV, this effect is 0.2%. Hence, we included the term $e^{-E_{\pi}(L_{t}-t)}$ in the denominator of Eq. (10) for the extrapolation of $R(t,\tau,z;P_{z})$ for $P_{z}=0$, and we also checked that the effect of the periodicity of lattice was indeed negligible for any of the non-zero $P_{z}$ we used.

For the case of $P_{z}=0$ displayed in the top-most panels of Fig. 22, we have shown the data in two ways to make the fits and the extrapolated value easier to understand. The unfilled symbols are the data for $R(t,\tau;z,P_{z},\gamma_{t})$ defined as the ratio of $C_{3{\rm pt}}(t,\tau;z,P_{z})$ to $C_{2{\rm pt}}(t;P_{z})$, and the solid curves are the fits including the $e^{-E_{\pi}(L_{t}-t)}$ term in the denominator of Eq. (10). While the fits describe the data well, the trend in the data with increasing $t$ can be seen to be away from the extrapolated value. To make the reason clearer, we have shown the modified ratio of $C_{3{\rm pt}}(t,\tau;z,P_{z})$ to the two-point function without the wrap-around term, $C_{2{\rm pt}}(t;P_{z})-A_{0}e^{-E_{\pi}(L_{t}-t)}$, as the filled symbols. The dashed curves are now the fits using just Eq. (10). The values of the amplitude $A_{0}$ and the energy $E_{\pi}$ were obtained by the two-state fit as described in the main text. Now, the trend with increasing $t$ is clearer.