# Lattice QCD Determination of the Bjorken-$x$ Dependence of Parton
Distribution Functions at Next-to-next-to-leading Order

Xiang Gao<EMAIL_ADDRESS>Key Laboratory of Quark & Lepton Physics (MOE) and
Institute of Particle Physics, Central China Normal University, Wuhan 430079,
China  Andrew D. Hanlon  Physics Department, Brookhaven National Laboratory,
Bldg. 510A, Upton, New York 11973, USA  Swagato Mukherjee  Physics Department,
Brookhaven National Laboratory, Bldg. 510A, Upton, New York 11973, USA  Peter
Petreczky  Physics Department, Brookhaven National Laboratory, Bldg. 510A,
Upton, New York 11973, USA  Philipp Scior  Physics Department, Brookhaven
National Laboratory, Bldg. 510A, Upton, New York 11973, USA  Sergey Syritsyn
RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, New York
11973  Department of Physics and Astronomy, Stony Brook University, Stony
Brook, New York 11790  Yong Zhao<EMAIL_ADDRESS>Physics Division, Argonne
National Laboratory, Lemont, IL 60439, USA

###### Abstract

We report the first lattice QCD calculation of pion valence quark distribution
with next-to-next-to-leading order perturbative matching correction, which is
done using two fine lattices with spacings $a=0.04$ fm and $0.06$ fm and
valence pion mass $m_{\pi}=300$ MeV, at boost momentum as large as $2.42$ GeV.
As a crucial step to control the systematics, we renormalize the pion valence
quasi distribution in the recently proposed hybrid scheme, which features a
Wilson-line mass subtraction at large distances in coordinate space, and
develop a procedure to match it to the ${\overline{\mathrm{MS}}}$ scheme. We
demonstrate that the renormalization and the perturbative matching in
Bjorken-$x$ space yield a reliable determination of the valence quark
distribution for $0.03\lesssim x\lesssim 0.80$ with 5–20% uncertainties.

Understanding the hadron inner structure remains one of the top fundamental
questions in nuclear and particle physics. As the lightest hadrons in nature,
pions are the Nambu-Goldstone bosons of quantum chromodynamics (QCD), and
their quark and gluon structures can help understand the origins of hadron
mass and dynamical chiral symmetry breaking. The parton distribution functions
(PDFs), which describe 1D momentum densities of quarks and gluons in a hadron,
are the simplest and most important quantities that have been extensively
studied from global high-energy scattering experiments and will be probed at
unprecedented precision at the future Electron-Ion Collider [1, 2]. Besides
the experimental efforts, the first-principles calculations of PDFs using
lattice QCD are also expected to provide useful predictions.

Computation of the PDFs on a Euclidean lattice has been extremely difficult
because they are defined from light-cone correlations with real-time
dependence in Minkowski space. For a long time, only the lowest moments of the
PDFs were calculable as they are matrix elements of local gauge-invariant
operators. For reviews see Refs. [3, 4]. Less than a decade ago, a
breakthrough was made by large-momentum effective theory (LaMET) [5, 6, 7],
which starts from a Euclidean “quasi-PDF” (qPDF) in a boosted hadron and
obtains the PDF through a large-momentum expansion and perturbative matching
of the qPDF in Bjorken-$x$ (longitudinal momentum fraction) space. Over the
years, LaMET has led to much progress in the calculation of PDFs and other
parton physics [7, 4], which reinvigorated the field as other proposals [8, 9,
10, 11, 12, 13] are also being studied and implemented.

Despite substantial progress, lattice calculation of the PDF $x$-dependence
has yet to achieve essential control of the systematic uncertainties [14]. In
the LaMET approach, lattice renormalization is one of the most important
sources of error. The nonlocal quark bilinear operator
$O_{\Gamma}(z)\equiv\bar{\psi}(z)\Gamma W(z,0)\psi(0)$, where $\Gamma$ is a
Dirac matrix and $z^{\mu}=(0,0,0,z)$, which defines the qPDF, suffers from a
linear power divergence in the Wilson line $W(z,0)$ that must be subtracted
before taking the continuum limit. The most popular methods so far are the
regularization independent momentum subtraction scheme [15, 16, 17, 18] and
other ratio schemes [19, 20, 21, 22], which use the matrix element of
$O_{\Gamma}(z)$ in an off-shell quark [15, 16, 17, 18], a static/boosted
hadron [19, 22] or the vacuum state [20, 21] as the renormalization factor. At
small $z$ the matrix elements in these schemes satisfy a factorization
relation to the light-cone correlation [23, 24, 13, 25]. However, at large $z$
they introduce nonperturbative effects [26] that propagate to the qPDF via
Fourier transform (FT) of the matrix elements, which contaminates the LaMET
matching in $x$-space. To overcome this limitation, the hybrid scheme [27] was
proposed to subtract the linear divergence at large $z$ and match the result
to the ${\overline{\mathrm{MS}}}$ scheme, thus preserving the LaMET matching
after FT. To date, the hybrid scheme has not been used in calculating the
PDFs, except for a recent work on meson distribution amplitudes [28]. Apart
from renormalization, the accuracy of perturbative matching also controls the
precision of the calculation. In all the existing lattice calculations, the
matching was done at only next-to-leading order (NLO), and it is not until
recently that the next-to-next-to-leading order (NNLO) matching was derived
for the non-singlet quark qPDF in the ${\overline{\mathrm{MS}}}$ scheme [29,
21].

In this Letter we present a state-of-the-art calculation of pion valence quark
PDF using high-statistics, superfine-spacing, and large-momentum lattice data
[26], with an adapted hybrid-scheme renormalization and the first-time
implementation of NNLO matching. The pion valence PDF has been extracted from
global fits [30, 31, 32, 33] and studied in lattice QCD [34, 35, 36, 37, 38,
39, 26, 40], with both at NLO accuracy. In this work, we subtract the linear
divergence in $O_{\Gamma}(z)$ with sub-percent precision, and develop a
procedure to match the lattice subtraction scheme to
${\overline{\mathrm{MS}}}$, a crucial step in the hybrid scheme to reduce the
power corrections [27]. We derive the NNLO hybrid-scheme matching and apply it
to the qPDF, showing good perturbative convergence and reduced scale-variation
uncertainty compared to NLO matching. Finally, we demonstrate that our
analysis yields a reliable determination of the PDF for $0.03\lesssim
x\lesssim 0.80$ with 5–20% uncertainties.

Our lattice data was produced using gauge ensembles in 2+1 flavor QCD
generated by the HotQCD collaboration [41] with Highly Improved Staggered
Quarks [42], including two lattice spacings $a=0.04$ and $0.06$ fm, and
volumes $L_{s}^{3}\times L_{t}=64^{4}$ and $48^{3}\times 64$, respectively. We
use tadpole-improved clover Wilson valence fermions on the hypercubic (HYP)
smeared [43] gauge background, with a valence pion mass $m_{\pi}=300$ MeV.
Furthermore, the Wilson line in $O_{\Gamma}(z)$ is constructed from HYP-
smeared gauge links. We use pion momenta $P^{z}=(2\pi n_{z})/(L_{s}a)$ with
$0\leq n_{z}\leq 5$, resulting in $P^{z}$ as large as $2.42$ GeV.

The qPDF $\tilde{f}_{v}(x,P^{z},\mu)$ is defined in a boosted pion state
$|P\rangle$ with four-momentum $P^{\mu}=(P^{t},0,0,P^{z})$:

$\displaystyle\tilde{f}_{v}(x,P^{z},\mu)$ $\displaystyle=\int{dz\over
2\pi}e^{ixP^{z}z}\ \tilde{h}(z,P^{z},\mu)\,,$ (1)

where $\tilde{h}(z,P^{z},\mu)\equiv\langle
P|O_{\gamma^{t}}(z)|P\rangle/(2P^{t})$, and $\mu$ is the
${\overline{\mathrm{MS}}}$ scale. The operator $O_{\Gamma}(z)$ can be
renormalized under lattice regularization as [44, 45, 46]

$\displaystyle O^{B}_{\Gamma}(z,a)=e^{-\delta
m(a)|z|}Z_{O}(a)O^{R}_{\Gamma}(z)\,,$ (2)

where “$B$” and “$R$” denote bare and renormalized quantities. The factor
$Z_{O}(a)$ includes all the logarithmic ultraviolet (UV) divergences which are
independent of $z$, while the Wilson-line mass correction $\delta m(a)$
includes the linear UV divergence $\propto 1/a$ and can be expressed as

$\displaystyle\delta m(a)={m_{-1}(a)\over a}+m_{0}\,,$ (3)

where $m_{-1}(a)$ is a series in the strong coupling $\alpha_{s}(1/a)$, and
$m_{0}$ is an ${\cal O}(\Lambda_{\rm QCD})$ constant originating from the
renormalon ambiguity in $m_{-1}(a)$ [47].

The hybrid scheme is implemented as follows: For $0\leq z\leq z_{S}$ with
$a\ll z_{S}\ll 1/\Lambda_{\rm QCD}$, we form the ratio
$\tilde{h}(z,P^{z},a)/\tilde{h}(z,0,a)$ to cancel the UV divergences and the
cutoff effects from $z\sim a$ [19]; at $z>z_{S}$ we subtract $\delta m(a)$ and
determine $Z_{O}(a)$ by imposing a continuity condition of the renormalized
matrix elements at $z=z_{S}$. There are different ways to calculate $\delta
m(a)$ [27, 48, 46, 49, 50, 51]. We determine $\delta m(a)$ from the
combination of the static quark-antiquark potential, $V^{\rm lat}(r)$ [41,
52], and the free energy of a static quark at non-zero temperature [53, 54,
55], with the following normalization scheme,

$\displaystyle V^{\rm lat}(a,r=r_{0})+2\delta m(a)=0.95/r_{0}\,,$ (4)

where $r_{0}=0.469$ fm is the Sommer scale for 2+1 flavor QCD [41], and the
constant 0.95 defines the scheme. The linear divergence $m_{-1}(a)/a$ does not
depend on the scheme, while $m_{0}$ does. The results are $a\delta
m=0.1586(8)$ and $0.1508(12)$ for $a=0.06$ and $0.04$ fm, respectively.

Since $m_{0}$ is scheme dependent, a factor of $e^{\overline{m}_{0}|z|}$ with
$\overline{m}_{0}\sim{\cal O}(\Lambda_{\rm QCD})$ is needed to match the
lattice scheme to ${\overline{\mathrm{MS}}}$, otherwise the LaMET expansion of
the qPDF will include a power correction $\propto\bar{m}_{0}/P^{z}$ [27],
which slows down convergence to the PDF as $P^{z}$ grows. It was proposed that
$\overline{m}_{0}$ can be obtained by comparing the subtracted matrix elements
of $O_{\Gamma}(z)$ [51] or $W(z,0)$ [49] with their ${\overline{\mathrm{MS}}}$
operator product expansion (OPE), whose accuracy requires $z\lesssim 0.2$ fm
[27]. But due to discretization effects, the window of $z$ that can be used is
actually narrow.

Our new procedure for the hybrid scheme is distinct by the determination of
$\overline{m}_{0}$. In order to use larger $z$, we construct the following
ratio and compare it to a form motivated by the OPE of $\tilde{h}(z,0,\mu)$,

$\displaystyle\lim_{a\to 0}e^{\delta
m(a)(z-z_{0})}{\tilde{h}(z,0,a)\over\tilde{h}(z_{0},0,a)}$
$\displaystyle\\!=\\!e^{-{\bar{m}_{0}}(z-z_{0})}{C_{0}(\mu^{2}z^{2})+\Lambda
z^{2}\over C_{0}(\mu^{2}z^{2}_{0})+\Lambda z^{2}_{0}}\,,$ (5)

where $z,z_{0}\gg a$, and the parameter $\Lambda\sim{\cal O}(\Lambda_{\rm
QCD}^{2})$. The Wilson coefficient $C_{0}$ is known to NNLO [25, 29, 21], and
$\bar{m}_{0}$ and $\Lambda z^{2}$ originate from the leading UV and infrared
renormalons in $C_{0}$ [20]. According to Eq. (2), the l.h.s. of Eq. (5) must
have a continuum limit if $\delta m(a)$ includes all the linear divergences,
which is renormalization group (RG) invariant. We choose $z\geq z_{0}=0.24$ fm
and find agreement between the $a=0.04$ fm and $a=0.06$ fm ratios at sub-
percent level up to $z\sim 1$ fm (see App. A). Then we extrapolate the lattice
ratios to the continuum with $a^{2}$-dependence [26], and fit the result to
the r.h.s of Eq. (5). For $z_{0}\leq z\leq 0.4$ fm, we obtain decent plateaus
and $\chi^{2}$ values for both $\bar{m}_{0}$ and $\Lambda$ with the NNLO
$C_{0}$. By definition $\bar{m}_{0}$ cancels the lattice scheme dependence of
$\delta m(a)$, as changing the scheme only shifts $\delta m(a)$ by a constant,
but $\bar{m}_{0}$ will inherit the ${\cal O}(\Lambda_{\rm QCD})$ ambiguity in
the ${\overline{\mathrm{MS}}}$ scheme. Since $C_{0}$ is at fixed order, both
$\bar{m}_{0}$ and $\Lambda$ depend on $\mu$, which we vary to estimate the
related uncertainty in the final result. At $\mu=2.0$ GeV,
$\bar{m}_{0}=0.151(1)$ GeV and $\Lambda=0.041(6)$ GeV2, so the power
correction is not negligible. Therefore, we modify the hybrid scheme by
correcting the $\Lambda z^{2}$ term in $\tilde{h}(z,0,\mu)$ at short $z$ as

$\displaystyle\\!\tilde{h}(z,z_{S},P^{z},\mu,a)\\!=\\!N{\tilde{h}(z,P^{z},a)\over\tilde{h}(z,0,a)}\\!{C_{0}(z^{2}\mu^{2})\\!+\\!\Lambda
z^{2}\over C_{0}(z^{2}\mu^{2})}\theta(z_{S}\\!-\\!z)$
$\displaystyle+\\!Ne^{\delta
m^{\prime}(z-z_{S})}{\tilde{h}(z,P^{z},a)\over\tilde{h}(z_{S},0,a)}{C_{0}(z^{2}_{S}\mu^{2})\\!+\\!\Lambda
z^{2}_{S}\over C_{0}(z^{2}_{S}\mu^{2})}\theta(z\\!-\\!z_{S})\,,\\!$ (6)

where $\delta m^{\prime}=\delta m+\bar{m}_{0}$, and
$N=\tilde{h}(0,0,a)/\tilde{h}(0,P^{z},a)$ normalizes
$\tilde{h}(z,z_{S},P^{z},\mu,a)$ to one at $z=0$. Since $C_{0}$ is at fixed
order, $\tilde{h}(z,z_{S},P^{z},\mu,a)$ depends on $\mu$ despite the fact that
it should be RG invariant. Such a renormalization is performed through
bootstrap loops so that the correlation between different $P^{z}$ and $z$ is
taken care of.

The hybrid-scheme matrix elements are shown in Fig. 1. At small $z$,
$\tilde{h}(z,P^{z})$ is dominated by the leading-twist contribution. At large
$z$, the spacelike correlator for pion valence quarks will exhibit an
exponential decay $\propto e^{-m_{\rm eff}|z|}$ where $m_{\rm eff}$ is an
effective mass related to the system [56]. When plotted as a function of
$\lambda=zP^{z}$, $\tilde{h}(\lambda,P^{z})$ should scale in $P^{z}$ at small
$\lambda$, with slight violation due to QCD evolution. Its exponential decay
will emerge at a larger $\lambda$ with greater $P^{z}$ and with decay rate
$m_{\rm eff}/P^{z}$. In the $P^{z}\to\infty$ limit, the exponential decay
vanishes at finite $\lambda$ ($z\to 0$), and only the leading-twist
contribution remains, which almost scales in $P^{z}$ and features a power-law
decay at large $\lambda$ that corresponds to small-$x$ PDF [27]. This picture
is consistent with Fig. 1.

Figure 1: Renormalized matrix elements in the hybrid scheme.

The next step is a FT. We truncate the matrix elements at $z_{L}$ or
$\lambda_{L}=z_{L}P^{z}$ where $\tilde{h}(\lambda_{L})\sim 0$, and extrapolate
to $\infty$ to remove the unphysical oscillations from a truncated FT [27].
The extrapolation form is $Ae^{-m_{\rm eff}|z|}/|\lambda|^{d}$, where $A$,
$m_{\rm eff}$ and $d$ are the parameters. Since $m_{\rm eff}$ is independent
of $P^{z}$, by fitting to the $P^{z}=0$ matrix elements we find that it is
around $0.1$ GeV, which is not far from the phenomenological estimate of
0.2–0.5 GeV in HQET [57]. Therefore, we impose $m_{\rm eff}>0.1$ GeV, as well
as $A>0$ and $d>0$, to ensure a convergent FT on each bootstrap sample. Since
the FT converges fast with the exponential decay, the extrapolation mainly
affects the small-$x$ region apart from removing the unphysical oscillations.
To verify this we vary $z_{L}$, which turns out to have little impact, and use
different $m_{\rm eff}$ bounds and extrapolation forms, which lead to
consistent qPDFs down to $x\sim 0.05$. (See App. B).

Then, we match the qPDF $\tilde{f}_{v}(x,\lambda_{S},P^{z},\mu)$ to the
${\overline{\mathrm{MS}}}$ PDF $f_{v}(x,\mu)$ through LaMET [58, 59, 25, 27]:

$\displaystyle f_{v}(x,\mu)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dy}{|y|}\
C^{-1}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)\tilde{f}_{v}(y,\lambda_{S},P^{z},\mu)$
$\displaystyle\qquad+{\cal
O}\Big{(}\frac{\Lambda_{\text{QCD}}^{2}}{(xP^{z})^{2}},\frac{\Lambda_{\text{QCD}}^{2}}{((1-x)P^{z})^{2}}\Big{)}\,,$
(7)

where $\lambda_{S}=z_{S}P^{z}$, $z_{S}=0.24$ fm, and the power corrections are
controlled by the parton and spectator momenta $xP^{z}$ and $(1-x)P^{z}$ [27].
Here $C^{-1}$ is the inverse of the hybrid-scheme matching coefficient $C$,
which we derive at NNLO [60] by conversion from the ${\overline{\mathrm{MS}}}$
result [29, 21]. Based on Eq. (Lattice QCD Determination of the Bjorken-$x$
Dependence of Parton Distribution Functions at Next-to-next-to-leading Order),
we can directly calculate the PDF with $P^{z}$-controlled power corrections
for $x\in[x_{\rm min},x_{\rm max}]$.

Figure 2: Comparison of PDFs obtained from the qPDF with NLO and NNLO matching
corrections.

In Fig. 2 we show the results of perturbative matching. The matching drives
the qPDF to smaller $x$ and reduces the statistical errors at moderate $x$,
because matching effectively relates the qPDF from finite $P^{z}$ to infinity,
and the qPDF evolves to smaller $x$ as $P^{z}$ increases. The NNLO correction
is generally smaller than the NLO correction, which indicates good
perturbative convergence, a crucial criterion for precision calculation.
Besides, by varying $\mu$ and evolving the matched results to the same $\mu$,
we find that the scale-variation uncertainty is reduced at NNLO, which is
further evidence of improved precision. The matching correction diverges as
$x\to 0$, implying that resummation of small-$x$ logarithms is needed. A
resummation is also necessary as $x\to 1$ [40], but these resummations are not
needed for moderate $x$.

We compare the PDFs obtained at different $P^{z}$ with NNLO matching in Fig.
3. At moderate $x$, the $P^{z}$-dependence is remarkably reduced, and the
results appear to converge for $P^{z}\geq 1.45$ GeV, which strongly indicates
the effectiveness of LaMET matching. At $x\gtrsim 1$, each PDF curve has a
small non-vanishing tail due to the power corrections in Eq. (Lattice QCD
Determination of the Bjorken-$x$ Dependence of Parton Distribution Functions
at Next-to-next-to-leading Order), but they decrease with larger $P^{z}$ (see
also App. C.3). To estimate the size of the power corrections, we fit the PDFs
obtained at $a=0.04$ fm, $P^{z}=\\{1.45,1.94,2.42\\}$ GeV and $a=0.06$ fm,
$P^{z}=\\{1.72,2.15\\}$ GeV to the ansatz $f_{v}(x)+\alpha(x)/P_{z}^{2}$ for
each fixed $x$, where we ignore the $a$-dependence as it has been found that
the matrix elements have ${\cal O}(a^{2}P_{z}^{2})$ effects that are less than
1% [26]. Since this fit is mainly affected by the data sets at lower $P^{z}$
with smaller statistical errors, which have larger power corrections, we use
the result at $P^{z}=2.42$ GeV instead of the fitted $f_{v}(x)$ as our final
prediction. The power correction at $P^{z}=2.42$ GeV is estimated to be
$\alpha(x)/[P_{z}^{2}f_{v}(x)]<0.10$ for $0.01<x<0.80$. It is surprising that
the results are insensitive to $P^{z}$ for $x$ as small as $0.01$, nor do they
show dependence on the extrapolation form in the FT as we have checked. This
can be explained by that, under matching, the qPDF contributes to the PDF at
larger $x$ which has less dependence on $P^{z}$ or the extrapolation.
Nevertheless, it must be pointed out that the smallness here is only relative,
as $\alpha(x)/P_{z}^{2}$ still diverges as $x\to 0$.

Figure 3: The PDFs obtained from the qPDFs with NNLO matching at different
$P^{z}$.

Our final prediction for the pion valence quark PDF (BNL-ANL21) is shown in
Fig. 4, which is obtained from the qPDF at $a=0.04$ fm, $z_{S}=0.24$ fm,
$z_{L}=0.92$ fm, $\mu=2.0$ GeV and $P^{z}=2.42$ GeV with exponential
extrapolation and NNLO matching. The red band represents the statistical
error, and the light purple band includes the error from scale variations,
which is obtained by repeating the same analysis for $\mu=1.4$ GeV and $2.8$
GeV and evolving the PDFs to $\mu=2.0$ GeV with the NLO DGLAP kernel. Since
the hybrid-scheme parameter $\bar{m}_{0}$ depends on $\mu$, the small scale
variation in the final result shows that the renormalization uncertainty is
well under control. We require that the ${\cal O}(\alpha_{s}^{3})$ matching
correction at $\mu=2.0$ GeV be smaller than 5%, which propagates geometrically
to $<37\%$ at NLO and $<14\%$ at NNLO, thus excluding $x<0.03$ and $x>0.88$. A
list of the above uncertainties at selected $x$ is shown in Table 1. See also
App. D. We neglect the FT uncertainty as it is extremely small. As for
$m_{\pi}$ dependence, our associated calculation of the second PDF moment at
$m_{\pi}=140$ MeV [61] shows consistency within 5% statistical uncertainty,
which will be validated by a direct comparison in the future. Previous studies
[62, 63] also suggest that the finite volume correction is less than 1% for
our lattice setup. At last, by limiting the estimated power corrections to be
less than 10%, we determine the PDF at $0.03\lesssim x\lesssim 0.80$ with
5–20% uncertainties. Our result is in great agreement with the recent global
fits by xFitter [31] and JAM21nlo [32] for $0.2<x<0.6$, but deviates from the
earlier GRVPI1 [30] and ASV [33] fits. When compared to a previous analysis of
the same lattice data (BNL20) [26], which used a short-distance factorization
of the matrix elements at NLO, and a parameterization of the PDF, our new
result has shifted central values and considerably reduced uncertainties at
moderate $x$, but still agrees within errors. With finite $P^{z}$ and
statistics, lattice QCD can only make predictions for $x\in[x_{\rm min},x_{\rm
max}]$. The PDF parameterization correlates the information at all
$x\in[0,1]$, so the larger uncertainties at moderate $x$ in BNL20 could be
propagated from the uncontrolled errors in the end-point regions. Besides,
there is no practical estimate of the model uncertainty in the
parameterization. Therefore, the LaMET calculation for $x\in[x_{\rm
min},x_{\rm max}]$ is more reliable as it does the power expansion and
matching directly in $x$-space.

Figure 4: Comparison of our prediction of $f_{v}(x)$, BNL-ANL21, to global fits and BNL20. The shaded regions $x<0.03$ and $x>0.8$ are excluded by requiring that estimates of ${\cal O}(\alpha_{s}^{3})$ and power corrections be smaller than 5% and 10%, respectively. $x$ | Statistical | Scale | ${\cal O}(\alpha_{s}^{3})$ | Power corrections | ${\cal O}(a^{2}P_{z}^{2})$
---|---|---|---|---|---
0.03 | 0.10 | 0.04 | $<0.05$ | $<0.01$ | $<0.01$
0.40 | 0.07 | $<0.01$ | $<0.05$ | 0.04 | $<0.01$
0.80 | 0.15 | 0.03 | $<0.05$ | 0.10 | $<0.01$
Table 1: Statistical and systematic uncertainties at given $x$.

In summary, we have performed a state-of-the-art lattice QCD calculation of
the $x$-dependence of pion valence quark PDF, where we developed a procedure
to renormalize the qPDF in the hybrid scheme and match it to the
${\overline{\mathrm{MS}}}$ PDF at NNLO. The final results show reduced
perturbation theory uncertainty and converge at moderate $x$ with pion momenta
greater than $1.45$ GeV, which allows us to reliably estimate the systematic
errors. This calculation can be improved with physical pion mass, continuum
extrapolation, and higher statistics for the matrix elements at long distances
and at larger boost momenta.

Our renormalization procedure can also be incorporated into the lattice
calculations of gluon PDFs, distribution amplitudes, generalized parton
distributions and transverse momentum distributions. With the systematics
under control, we can expect lattice QCD to provide reliable predictions for
these quantities in the future.

###### Acknowledgements.

We thank Vladimir Braun, Xiangdong Ji, Nikhil Karthik, Yizhuang Liu, Antonio
Pineda, Yushan Su and Jianhui Zhang for valuable communications. This material
is based upon work supported by: (i) The U.S. Department of Energy, Office of
Science, Office of Nuclear Physics through Contract No. DE-SC0012704 and No.
DE-AC02-06CH11357; (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 U.S. Department of Energy, Office of Science, Office of
Nuclear Physics, within the framework of the TMD Topical Collaboration. (iv)
XG is partially supported by the NSFC under the grant number 11775096 and the
Guangdong Major Project of Basic and Applied Basic Research No.
2020B0301030008. (v) SS is supported by the National Science Foundation under
CAREER Award PHY-1847893 and by the RHIC Physics Fellow Program of the RIKEN
BNL Research Center.. (vi) This research used awards of computer time provided
by the INCITE and ALCC programs at Oak Ridge Leadership Computing Facility, a
DOE Office of Science User Facility operated under Contract No. DE-
AC05-00OR22725. (vii) Computations for this work were carried out in part on
facilities of the USQCD Collaboration, which are funded by the Office of
Science of the U.S. Department of Energy. viii) YZ is partially supported by
an LDRD initiative at Argonne National Laboratory under Project No. 2020-0020.

##  Appendix A Hybrid scheme renormalization

###  A.1 Definition of scheme

As has been described in the main text, the hybrid scheme renormalization
includes two parts:

  * •

For $z\leq z_{S}$, we form the ratio of bare matrix elements [19],

$\displaystyle\frac{\tilde{h}(z,P^{z},a)}{\tilde{h}(z,0,a)}\,,$ (8)

which has a well-defined continuum limit and is renormalization group (RG)
invariant.

  * •

For $z\geq z_{S}$, the renormalized matrix element is

$\displaystyle e^{\delta
m(a)|z-z_{S}|}{\tilde{h}(z,P^{z},a)\over\tilde{h}(z_{S},0,a)}\,,$ (9)

which is equal to the ratio in Eq. (8) at $z=z_{S}$. To determine $\delta
m(a)$ we use the additive renormalization constant, $c_{Q}(a)=\delta m(a)$,
which is obtained in Ref. [54] from the analysis of the free energy of a
static quark, $F_{Q}(T)$, at non-zero temperature $T$ with the normalization
condition in Eq. (4). Recently $F_{Q}$ has been calculated using one step of
HYP smearing [55], and it was found that HYP smearing does not affect the
temperature dependence of $F_{Q}(T)$, but only shifts it by an additive
constant. Therefore, we have $F_{Q}^{B,1}(T)+\delta
m(a)=F_{Q}^{B,0}(T)+c_{Q}(a)$ with superscripts 0 and 1 referring to the
number of HYP smearing steps in the bare free energy of the static quark.
Using the lattice results for $F_{Q}^{B,0}(T)$ and $F_{Q}^{B,1}(T)$ obtained
on $N_{\tau}=12$ lattices and temperatures corresponding to $a=0.04$ fm and
$a=0.06$ fm (where cutoff effects can be neglected), as well as the values of
$c_{Q}$ from Table X of Ref. [54] for $\beta=7.825$ ($a=0.04$ fm) and
$\beta=7.373$ ($a=0.06$ fm), we obtain $\delta m(a)$. The results are $a\delta
m(a=0.06{\rm\ fm})=0.1586(8)$ and $a\delta m(a=0.04{\rm\ fm})=0.1508(12)$.

First of all, to test how well the subtraction of $\delta m(a)$ can remove the
linear divergences in $\tilde{h}(z,P^{z},a)$, we construct the ratio in Eq.
(5),

$\displaystyle\tilde{R}(z,z_{0},a)$ $\displaystyle\equiv e^{\delta
m(a)(z-z_{0})}{\tilde{h}(z,0,a)\over\tilde{h}(z_{0},0,a)}\,,$ (10)

where $z_{0}=0.24$ fm for both lattice spacings. According to Eq. (2), the
renormalization factor $Z_{O}(a)$ cancels out in the ratio. Therefore, if
$\delta m(a)$ includes all the linear divergences, then $\tilde{R}(z,z_{0},a)$
should have a well-defined continuum limit.

Figure 5: Upper panel: ratios of bare lattice matrix elements without the
Wilson-line mass subtraction. Lower panel: the ratio in Eq. (10) with Wilson-
line mass subtraction. The red and blue points are for $a=0.04$ fm and $0.06$
fm. The red and blue bands are interpolations of the points, and the gray band
is the continuum extrapolation of them with $a^{2}$-dependence.

Our lattice results for the above ratio with $z_{0}=0.24$ fm is shown in Fig.
5. As one can see, the differences between the ratios at $a=0.04$ fm and
$0.06$ fm are at sub-percent level, which clearly shows that the linear
divergences have been sufficiently subtracted by $\delta m(a)$. Therefore, the
ratio in Eq. (10) has a continuum limit

$\displaystyle\lim_{a\to 0}\tilde{R}(z,z_{0},a)$
$\displaystyle=\tilde{R}(z,z_{0})\,,$ (11)

which is RG invariant.

Our next step is to match the lattice subtraction scheme to
${\overline{\mathrm{MS}}}$. When $z,z_{0}\ll\Lambda_{\rm QCD}^{-1}$, the
${\overline{\mathrm{MS}}}$ matrix element
$\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$ has an OPE that goes as

$\displaystyle\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$
$\displaystyle=e^{-m^{{\overline{\mathrm{MS}}}}_{0}|z|}\left[C_{0}(z^{2}\mu^{2})\right.$
$\displaystyle\left.+z^{2}C_{2}(z^{2}\mu^{2})\langle P|O_{\rm
tw4}(\mu)|P\rangle+\ldots\right]\,,$ (12)

where $m^{{\overline{\mathrm{MS}}}}_{0}$ is the ${\cal O}(\Lambda_{\rm QCD})$
renormalon ambiguity from the Wilson line self-energy renormalization [20],
$O_{\rm tw4}(\mu)$ is a twist-four operator (for example,
$\bar{\psi}D^{2}\psi$ or $g\bar{\psi}\sigma_{\mu\nu}F^{\mu\nu}\psi$), $C_{0}$
and $C_{2}$ are perturbative coefficient functions, and “$\ldots$” denotes
contributions at higher twists. Since $P^{z}=0$, $C_{0}$ is the only Wilson
coefficient that contributes at leading-twist. The leading-twist contribution
is proportional to $\langle P|\bar{\psi}\gamma^{t}\psi|P\rangle/(2P^{t})$
which is trivially one due to vector current conservation. Since
$\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$ is multiplicatively
renormalizable, both $C_{0}(z^{2}\mu^{2})$ and $C_{2}(z^{2}\mu^{2})\langle
P|O_{\rm tw4}(\mu)|P\rangle$ must satisfy RG equations with the same anomalous
dimension, which is known to next-to-next-to-next-to-leading order (N3LO)
[64]. Due to the ambiguity in summing the perturbative series in
$C_{0}(z^{2}\mu^{2})$, there are $O(\Lambda_{\rm QCD}^{2n})$ IR renormalons in
the leading-twist contribution that should be cancelled by those from higher-
twist condensates, along with the ${\cal O}(\Lambda_{\rm QCD})$ UV renormalon
to be cancelled by $m^{{\overline{\mathrm{MS}}}}_{0}$ [20, 57]. Both the UV
and IR renormalon contributions cannot be well defined unless one specifies
how to sum the perturbative series in $C_{0}(z^{2}\mu^{2})$ to all orders,
which, however, is unknown as $C_{0}(z^{2}\mu^{2})$ has been calculated to
only NNLO so far [21].

Note that $m^{{\overline{\mathrm{MS}}}}_{0}$ is analogous to the mass
renormalization in heavy-quark effective theory (HQET) [57], which is of UV
origin and cannot be attributed to any short-distance condensate. Instead, it
appears as a residual mass term in the HQET Lagrangian and exists in
$\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$ at all $z$, i.e.,

$\displaystyle\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$
$\displaystyle=e^{-m^{{\overline{\mathrm{MS}}}}_{0}|z|}\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,0,\mu)\,,$
(13)

where $\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,0,\mu)$ at short distance
reduces to the OPE series in the brackets of Eq. (A.1).

The renormalons have been studied extensively for the Polyakov loop and
plaquette in lattice QCD [47, 65, 66, 67, 68]. In lattice perturbation theory,
one has to compute the perturbative series to very high orders of $\alpha_{s}$
in order to see the renormalon effects. Nevertheless, in the
${\overline{\mathrm{MS}}}$ scheme, the OPE with Wilson coefficient at a few
loop orders and the condensate term, turns out to be successful in describing
the static potential at short distance up to $\sim 0.25$ fm [69]. One
explanation is that $\alpha_{s}$ in the ${\overline{\mathrm{MS}}}$ scheme is
larger than that in lattice perturbation theory, so the renormalon effect
which is of ${\cal O}(\alpha_{s}^{n})$ with
$n\sim(2\pi)/(\beta_{0}\alpha_{s})$ becomes significant at lower orders. This
situation is similar to the OPE in QCD sum rules [70, 71, 72, 73, 74], which
works well in phenomenology. The reason behind such success is probably due to
a proper choice of the renormalization scale $\mu$ so that $\alpha_{s}(\mu)$
is small enough for the perturbative series to converge, while the
$\mu$-dependent effects in the condensate remain insignificant as they should
be of the same magnitude of highest order in the truncated perturbative series
[73, 74].

Therefore, we approximate Eq. (A.1) as

$\displaystyle\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$
$\displaystyle\approx
e^{-m^{{\overline{\mathrm{MS}}}}_{0}(\mu)|z|}\left[C_{0}^{\rm
FO}(z^{2}\mu^{2})+\Lambda(\mu)z^{2}\right]\,,$ (14)

where “FO” stands for fixed order, $\Lambda(\mu)$ is a parameter of ${\cal
O}(\Lambda_{\rm QCD}^{2})$, and we ignore the higher power corrections by
working at not too large $z$. The $\mu$ dependence of the parameters
$m^{{\overline{\mathrm{MS}}}}_{0}$ and $\Lambda$ is understandable because
this approximation is valid for a small window of $\mu$, and they also depend
on the perturbative orders in $C_{0}^{\rm FO}$ if the latter does not converge
fast. Note that the although the model in Eq. (14) is not guaranteed to
satisfy the RG equation for $\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$,
we argue that within the range of $\mu$ where it can describe the physical
results, the $\mu$-dependence in the power correction term, which is already
suppressed, is weak and can be ignored.

Based on the above approximation, we fit our lattice results of the ratio in
Eq. (10) to the following ansatz,

$\displaystyle\tilde{R}(z,z_{0})$
$\displaystyle=e^{-\bar{m}_{0}(\mu)(z-z_{0})}\frac{C_{0}^{\rm
FO}(z^{2}\mu^{2})+\Lambda(\mu)z^{2}}{C_{0}^{\rm
FO}(z_{0}^{2}\mu^{2})+\Lambda(\mu)z_{0}^{2}}\,,$ (15)

where the mass shift

$\displaystyle\bar{m}_{0}(\mu)$
$\displaystyle=-m_{0}+m^{{\overline{\mathrm{MS}}}}_{0}(\mu)\,,$ (16)

cancels the lattice scheme dependence of $m_{0}$ in Eq. (3) and introduces the
renormalon ambiguity of the ${\overline{\mathrm{MS}}}$ scheme. Effectively,
$\bar{m}_{0}$ matches the hybrid-scheme matrix elements at $z\geq z_{S}$ to
the ratio of $\tilde{h}_{0}^{{\overline{\mathrm{MS}}}}$ as

$\displaystyle\lim_{a\to 0}e^{(\delta
m(a)+\bar{m}_{0}(\mu))(z-z_{S})}{\tilde{h}(z,P^{z},a)\over\tilde{h}(z_{S},0,a)}$
$\displaystyle={\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,P^{z},\mu)\over\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z_{S},0,\mu)}\,.$
(17)

Moreover, since the ansatz in Eq. (15) can describe the short-distance matrix
elements well, we can correct the $\Lambda z^{2}$ term in
$\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,0,\mu)$ at $z\leq z_{S}$ as

$\displaystyle{\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,0,\mu)}\frac{C_{0}^{\rm
FO}(z^{2}\mu^{2})}{C_{0}^{\rm FO}(z^{2}\mu^{2})+\Lambda(\mu)z^{2}}\,,$ (18)

which is equivalent to replacing
${\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,0,\mu)}$ by the perturbative
$C_{0}$, as in Eq. (Lattice QCD Determination of the Bjorken-$x$ Dependence of
Parton Distribution Functions at Next-to-next-to-leading Order). Eventually,
the continuum limit of the matched matrix element in Eq. (Lattice QCD
Determination of the Bjorken-$x$ Dependence of Parton Distribution Functions
at Next-to-next-to-leading Order) is

$\displaystyle\tilde{h}(z,z_{S},P^{z},\mu)$
$\displaystyle={\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,P^{z},\mu)\over
C_{0}^{\rm FO}(z^{2}\mu^{2})}\theta(z_{S}-|z|)$
$\displaystyle\qquad+{\tilde{h}^{{\overline{\mathrm{MS}}}}_{0}(z,P^{z},\mu)\over
C_{0}^{\rm FO}(z^{2}_{S}\mu^{2})}\theta(|z|-z_{S})\,,$ (19)

which is different from ${\overline{\mathrm{MS}}}$ through a perturbative
matching for all $z$ as long as $z_{S}\ll\Lambda_{\rm QCD}^{-1}$. Therefore,
the qPDF defined as FT of $\tilde{h}(z,z_{S},P^{z})$ is still factorizable.

Note that $\bar{m}_{0}(\mu)$ introduces the ambiguity
$m^{{\overline{\mathrm{MS}}}}_{0}(\mu)$ to the matched matrix elements.
Nevertheless, we argue that $C^{\rm FO}_{0}(\mu^{2}z^{2})$ at NNLO is
different from a particular summation prescription by ${\cal
O}(\alpha_{s}^{3})$ contributions, which cannot be smaller than the ambiguity
in $m^{{\overline{\mathrm{MS}}}}_{0}(\mu)$ as the latter reflects the
uncertainty in summing divergent perturbative series at sufficiently high
orders. Therefore, we can attribute the renormalon ambiguity in
$\bar{m}_{0}(\mu)$ to higher loop-order effects, and estimate the latter by
varying $\mu$ by a factor of $\sqrt{2}$ and $1/\sqrt{2}$. The range of $\mu$
we vary from cannot be too large. If $\mu$ is too small, then
$\alpha_{s}(\mu)$ becomes too large; if $\mu$ is too large, then we need to
resum the large $\ln(z^{2}\mu^{2})$ in $C_{0}(z^{2}\mu^{2})$. In both cases
the perturbative series converges slowly. In our analysis, we scan $\mu$
within $[0.9,2.0]$ GeV for $C_{0}^{\rm NLO}$ and $[1.4,3.2]$ GeV for
$C_{0}^{\rm NNLO}$ to study the scale dependence and uncertainty from
renormalon ambiguity.

###  A.2 Fitting of $\bar{m}_{0}$ and $\Lambda(\mu)$

Currently, the Wilson coefficient $C_{0}(\mu^{2}z^{2})$ is known to NNLO [29,
21] and its anomalous dimension has been calculated at three-loop order [64],

$\displaystyle
C_{0}\big{(}\mu^{2}z^{2},\alpha_{s}(\mu)\big{)}=1+a_{s}\left(2L+\frac{10}{3}\right)$
$\displaystyle+\\!a_{s}^{2}\left[\frac{13}{2}L^{2}\\!+\\!\frac{1461\\!+\\!28\pi^{2}}{54}L\\!+\\!\frac{38127\\!-\\!824\pi^{2}\\!-\\!4032\zeta(3)}{648}\right]$
$\displaystyle+a_{s}^{3}\left[\frac{143}{6}L^{3}+\Big{(}\frac{6127}{36}+\frac{91\pi^{2}}{27}\Big{)}L^{2}\right.$
$\displaystyle\left.\qquad+\frac{690939+760\pi^{4}-8976\pi^{2}-94068\zeta(3)}{972}L+400\right]$
$\displaystyle+O(a_{s}^{4})\,,$ (20)

where $a_{s}=\alpha_{s}/(2\pi)$, $L=\ln(\mu^{2}z^{2}/b_{0}^{2})$, and
$b_{0}=2e^{-\gamma_{E}}$,. The factor $400$ in the last square bracket is a
simple guess by assuming that the constant part of the perturbative correction
grows as a geometric series in the order of $a_{s}$.

We also consider the RG improved (RGI) Wilson coefficient [40]

$\displaystyle C_{0}^{\rm RGI}\big{(}\mu^{2},z^{2}\big{)}$
$\displaystyle=C_{0}\big{(}1,\alpha_{s}(b_{0}/z)\big{)}$ (21)
$\displaystyle\quad\times\exp\Big{[}\int_{b_{0}/z}^{\mu}{d\alpha_{s}(\mu^{\prime})}\
{\gamma_{\cal
O}(\alpha(\mu^{\prime}))\over\beta(\alpha_{s}(\mu^{\prime}))}\Big{]}\,,$

where $\gamma_{\cal O}$ is the anomalous dimension of the operator
$O_{\Gamma}(z,\mu)$, and
$\beta(\alpha_{s}(\mu))=d\alpha_{s}(\mu)/d\ln\mu^{2}$. In this way, we can
first factor out the evolution factor in Eq. (14) as it must be satisfied by
the full matrix element $\tilde{h}^{{\overline{\mathrm{MS}}}}(z,0,\mu)$, and
therefore construct the ratio $\tilde{R}(z,z_{0})$ in an explicitly
$\mu$-independent way.

Figure 6: The fixed-order and RGI Wilson coefficients $C_{0}(z^{2}\mu^{2})$ up
to N3LO.

We compare $C_{0}$ and $C_{0}^{\rm RGI}$ at NLO, NNLO and N3LO at $\mu=2.0$
GeV in Fig. 6. The strong coupling constants at each perturbative order are
defined by the corresponding $\Lambda_{\rm QCD}^{{\overline{\mathrm{MS}}}}$
with one-, two- and three- loop $\beta$ functions and $n_{f}=3$, which are
fixed by matching to $\alpha_{s}(\mu=2\ {\rm GeV})=0.293$. The latter is
obtained from $\Lambda_{\rm QCD}^{{\overline{\mathrm{MS}}}}=332$ MeV with
five-loop $\beta$-function and $n_{f}=3$, as has been calculated using the
same lattice ensembles [75]. As one can see, at $z>0.2$ fm the RGI Wilson
coefficients start to deviate significantly from the fixed-order ones, which
is mainly due to the large value of $\alpha_{s}$ as in RGI Wilson coefficients
as we evolve from $\mu$ to $1/\boldsymbol{z}$. This indicates that at $z>0.2$
fm, the scale uncertainty in the perturbative series is significant due to the
enhancement of non-perturbative effects, and to use OPE we should work at very
short distances ($z<0.2$ fm). However, there will not be enough room for
varying $z$ to satisfy $z\gg a$ so that discretization effects are suppressed.
Therefore, in our analysis we loosen our requirement for very small $z$ by
only using the ansatz in Eq. (15) and not considering the RGI Wilson
coefficients.

(a)

(b)

Figure 7:  Effective mass $\bar{m}_{0}^{\rm eff}(z)$ (a) and its slope
$\bar{m}_{2}^{\rm eff}(z)$ (b) vs $z$.

In Fig. 7a, we plot an effective mass $\bar{m}_{0}^{\rm eff}(z)$ which is
defined as

$\displaystyle\bar{m}_{0}^{\rm eff}(z)(z-z_{0})$
$\displaystyle\equiv-\ln{\tilde{h}(z,0,a)\over\tilde{h}(z_{0},0,a)}+\ln{C_{0}^{\rm
NNLO}(z^{2}\mu^{2})\over C_{0}^{\rm NNLO}(z^{2}_{0}\mu^{2})}\,,$ (22)

where $\mu=2.0$ GeV. If the twist-four condensate is negligible, then we
should expect a plateau in $z$, but Fig. 7a shows that it has an almost
constant nonzero slope at $z$ from $0.24$ fm up to $1.0$ fm. In Fig. 7b we
plot its slope

$\displaystyle\bar{m}_{2}^{\rm eff}(z)$ $\displaystyle={\bar{m}_{0}^{\rm
eff}(z)-\bar{m}_{0}^{\rm eff}(z-a)\over a}\,,$ (23)

which is consistent with being constant for a wide range of $z$. This suggests
that there is considerable quadratic $z$-dependence from the twist-four
condensate, as inclued in the ansatz in Eq. (15).

(a)

(b)

(c)

(d)

Figure 8:  Results for $\bar{m}_{0}(\mu)$ (a) and $\Lambda(\mu)$ (b) fitted
from $\tilde{R}(z,z_{0})$ for $z_{0}<z<z_{\rm max}$ and $z_{0}=0.24$ fm, with
NLO and NNLO Wilson coefficients at various values of $\mu$.

Our results for $\bar{m}_{0}$ and $\Lambda$ fitted from $\tilde{R}(z,z_{0})$
for $z_{0}<z<z_{\rm max}$ with $z_{0}=0.24$ fm are shown in Fig. 8. As one can
see, the two parameters remain constant in $z_{\rm max}$ up to around $0.5$ fm
within a small window of $\mu$, which is different with the NLO and NNLO
Wilson coefficients. At larger $z$, the higher-twist and
$\alpha_{s}\ln(z^{2}\mu^{2})$ effects become significant, which can no longer
be described by the simple ansatz in Eq. (15). In this work, we use
$\tilde{R}(z,z_{0})$ at $0.24$ fm $<z<0.4$ fm to fit the parameters at all
$\mu$ as input for the hybrid scheme renormalization and matching. To estimate
the uncertainty from the choice of $\mu$, we will match the qPDFs obtained at
different $\mu$ to the corresponding PDFs, and then evolve the final results
to $\mu=2.0$ GeV for comparison.

##  Appendix B Fourier transform (FT)

The qPDF is defined as the FT of $\tilde{h}(z,z_{S},P^{z})$ or
$\tilde{h}(\lambda,\lambda_{S},P^{z})$,

$\displaystyle\tilde{f}(x,z_{S},P^{z})$ $\displaystyle=\int{d\lambda\over
2\pi}\ e^{ix\lambda}\tilde{h}(\lambda,\lambda_{S},P^{z})\,.$ (24)

Since $\tilde{h}(\lambda,\lambda_{S},P^{z})$ is perturbatively matched from
the ${\overline{\mathrm{MS}}}$ scheme, the factorization formula should still
be valid for the corresponding qPDF $\tilde{f}(x,z_{S},P^{z})$ [27].
Therefore, we should integrate over all $z$ in the FT to obtain the
$x$-dependence of the qPDF. However, due to finite lattice size effects,
worsening signal-to-noise ratio and other systematics at large $z$, we have to
truncate $\tilde{h}(z,z_{S},P^{z})$ at $z=z_{L}$ and extrapolate to
$z\to\infty$ to complete the FT. As a result, the small-$x$ ($x\lesssim
1/\lambda_{L}$) region is the most sensitive to the extrapolation model, and
the corresponding systematic uncertainty cannot be well controlled. On the
other hand, the reliability of the $x\gtrsim 1/\lambda_{L}$ region depends on
the premises that the $\tilde{h}(z)$ is small at $z=z_{L}$ and exhibits an
exponential decay when $z_{L}$ is large enough. The first condition is easy to
understand as a truncated FT will lead to an unphysical oscillation in the
$x$-space with amplitude proportional to $|\tilde{h}(z_{L})|$, while the
exponential decay guarantees that the FT converges fast and the qPDF at
$x\gtrsim 1/\lambda_{L}$ has very little dependence on the specific model used
in the extrapolation.

In this section, we first derive that the equal-time correlator in a hadron
state does exhibit an exponential decay at large distances, then we
demonstrate that including this constraint in the extrapolation will lead to a
reliable FT in the moderate-to-large $x$ region. Finally, we perform the
extrapolated FT on our lattice results.

###  B.1 Matrix elements at large $z$

To begin with, let us consider a current-current correlation in the vacuum,
$\langle\Omega|J_{5}(x)J_{5}(0)|\Omega\rangle$, where
$J_{5}=\bar{q}\gamma_{5}q$ and $x^{2}<0$. If we ignore the existence of zero
modes and only consider gapped vacuum excitations, then

$\displaystyle\langle\Omega|J_{5}(x)J_{5}(0)|\Omega\rangle$
$\displaystyle=\sum_{n}\int{d^{3}k_{n}\over(2\pi)^{3}2E_{k_{n}}}\langle\Omega|J_{5}(x)|n\rangle\langle
n|J_{5}(0)|\Omega\rangle$
$\displaystyle=\sum_{n}\int{d^{3}k_{n}\over(2\pi)^{3}2E_{k}}\langle\Omega|J_{5}(0)|n\rangle\langle
n|J_{5}(0)|\Omega\rangle e^{-ix\cdot k_{n}}$
$\displaystyle=\sum_{n}|Z_{n}|^{2}\int{d^{4}k_{n}\over(2\pi)^{4}}{e^{-ix\cdot
k_{n}}\over k^{2}_{n}-m_{n}^{2}+i0}$ $\displaystyle=-{i\over
4\pi^{2}}\sum_{n}|Z_{n}|^{2}{m_{n}\over\sqrt{-x^{2}}}K_{1}(m_{n}\sqrt{-x^{2}})\,.$
(25)

where $Z_{n}$ is the overlap between the operator $J_{5}(x)$ and intermediate
sate $|n\rangle$. Here $m_{n}$ is the mass of the intermediate state particle,
and $K_{n}$ is the modified Bessel function of the second kind. Then, since

$\displaystyle\lim_{|x|\to\infty}{m_{n}\over\sqrt{-x^{2}}}K_{1}(m_{n}\sqrt{-x^{2}})$
$\displaystyle=\sqrt{{\pi\over 2}}{\sqrt{m_{n}}\over|x|^{3\over
2}}e^{-m_{n}|x|}\,,$ (26)

The correlation function should, therefore, be dominated by the exponential
decay of the lowest-lying state that overlaps with $J_{5}(x)$.

When the external state is a static hadron, it has also been shown that the
spacelike correlations exhibit an exponential decay at large distance [56].

We are interested in equal-time quark bilinear correlators in a boosted hadron
state, which can be expressed in terms of the product of two “heavy-light”
currents [44, 46], where the “heavy quark” $h_{\hat{x}}$ is an auxiliary field
defined along the $\hat{x}$ direction, similar to that in HQET.

Let us choose the external state to be a pion. According to Lorentz
covariance, we can decompose the correlation as

$\displaystyle\langle\pi(p)|\bar{q}(x)\gamma^{\mu}h_{\hat{x}}(x)\bar{h}_{\hat{x}}(0)q(0)|\pi(p)\rangle$
$\displaystyle=p^{\mu}f_{p}(p\cdot x,x^{2})+x^{\mu}f_{x}(p\cdot x,x^{2})\,,$
(27)

where the scalar functions $f_{p,x}(p\cdot x,x^{2})$ are analytic functions of
$p\cdot x$ and $x^{2}$. We can select the index $\mu$ such that $x^{\mu}=0$.
For example, we can choose $\mu=z$ when $x^{\mu}=(t,0,0,0)$ or $\mu=t$ when
$x^{\mu}=(0,0,0,z)$. The HQET corresponds to the timelike case, as

$\displaystyle p^{z}f_{p}(p\cdot x,x^{2})$
$\displaystyle=\sum_{n}\int{d^{3}k_{n}\over(2\pi)^{3}2E_{k_{n}}}e^{-ix\cdot(k_{n}-m_{Q}v-p)}$
$\displaystyle\quad\times\langle\pi(p)|\bar{q}\Gamma h_{v}|n\rangle\langle
n|\bar{h}_{v}q|\pi(p)\rangle\,,$ (28)

where $h_{v}$ is the effective heavy-quark field moving with velocity
$v^{\mu}$ and related to the QCD heavy quark $Q$ by the projection

$\displaystyle h_{v}(x)$ $\displaystyle=e^{im_{Q}v\cdot x}{1+\not{v}\over
2}Q(x)\,.$ (29)

The lowest intermediate state $|H(v)\rangle$ is a heavy-light meson with mass
$m_{H}=m_{Q}+\bar{\Lambda}$ and momentum $k^{\mu}=m_{H}v^{\mu}$, where $m_{Q}$
is the heavy quark pole mass, and $\bar{\Lambda}$ can be interpreted as the
mass of the constituent light quark or binding energy. Both $\bar{\Lambda}$
and $m_{Q}$ have ${\cal O}(\Lambda_{\rm QCD})$ renormalon ambiguities which
cancel between each other. In the $\Lambda_{\rm QCD}/m_{Q}\to 0$ limit,
$\bar{\Lambda}$ should be independent of the heavy quark mass, but can depend
on the light quark mass.

The matrix element $\langle\pi(p)|\bar{q}\Gamma h_{v}|H(v)\rangle$ is given by
the transition form factors [76],

$\displaystyle\langle\pi(p)|\bar{q}\Gamma h_{v}|H(v)\rangle$
$\displaystyle=-\mathrm{Tr}\left\\{\gamma_{5}\left[f_{1}(v\cdot
p)+f_{2}(v\cdot p){\not{p}\over v\cdot p}\right]\Gamma{\cal M}(v)\right\\}\,,$
(30)

where the form factors $f_{1}$ and $f_{2}$ only depend on $v\cdot p$ in HQET,
and the projetion operator ${\cal M}(v)$ depends on the spin of the heavy-
light meson $H(v)$,

$\displaystyle{\cal M}(v)$ $\displaystyle={1+\not{v}\over
2}\Bigg{\\{}\begin{array}[]{cc}-\gamma_{5}\,,&{\rm for}\ J^{P}=0^{-}\,,\\\
\not{\epsilon}\,,&{\rm for}\ J^{P}=1^{-}\,,\end{array}$ (33)

with $\epsilon^{\mu}$ being the polarization vector for vector mesons.
Therefore,

$\displaystyle\langle\pi(p)|\bar{q}\gamma^{\mu}h_{v}|H(v)\rangle$
$\displaystyle=2f_{1}(v\cdot p)v^{\mu}+2f_{2}(v\cdot p)p^{\mu}\,,$ (34)
$\displaystyle\langle\pi(p)|\bar{q}h_{v}|H(v)\rangle$
$\displaystyle=2f_{1}(v\cdot p)+2f_{2}(v\cdot p)\,.$ (35)

Then, the correlation function becomes

$\displaystyle p^{z}f_{p}(p\cdot x,x^{2})\approx
4m_{Q}^{2}\sum_{n}\int{d^{3}\vec{v}\over(2\pi)^{3}2\sqrt{1+\vec{v}^{2}}}$
$\displaystyle\qquad\times
e^{-i\left(\bar{\Lambda}\sqrt{1+\vec{v}^{2}}-p^{0}\right)x^{0}}(f_{1}+f_{2})(f_{1}v^{z}+f_{2}p^{z})\,.$
(36)

Note that when $x^{0}\to\infty$, $x^{0}\bar{\Lambda}\sqrt{1+\vec{v}^{2}}\geq
x^{0}\bar{\Lambda}$ constitutes a large phase, so the integrand is quickly
oscillating and should be suppressed. To have a naive estimate, let us assume
$f_{1}$ and $f_{2}$ are constant in $v\cdot p$, and the remaining integral is
simply

$\displaystyle\int{d^{3}\vec{v}\over(2\pi)^{3}2\sqrt{1+\vec{v}^{2}}}e^{-i\left(\bar{\Lambda}\sqrt{1+\vec{v}^{2}}-p^{0}\right)x^{0}}$
$\displaystyle\qquad={1\over
4\pi^{2}}K_{1}\left(\bar{\Lambda}\sqrt{-x_{0}^{2}}\right){e^{ip^{0}x^{0}}\over\bar{\Lambda}\sqrt{-x_{0}^{2}}}$
$\displaystyle\qquad={1\over
4\pi^{2}}K_{1}\left(\bar{\Lambda}\sqrt{-x^{2}}\right){e^{ip\cdot
x}\over\bar{\Lambda}\sqrt{-x^{2}}}\,,$ (37)

where we first obtained the result for imaginary $x^{0}$ and then analytically
continued back to the real axis.

Then, using Lorentz invariance and analyticity, we can obtain the result for
$x^{2}<0$, which corresponds to the equal-time correlator that we calculate in
this work. At large separation, we have

$\displaystyle\lim_{|x|\to\infty}f_{p}(p\cdot x,x^{2})$ $\displaystyle\propto
m_{Q}^{2}{e^{-\bar{\Lambda}|x|}\over(\bar{\Lambda}|x|)^{3\over 2}}e^{ip\cdot
x}\,,$ (38)

which also exhibits an exponential behavior with decay constant
$\bar{\Lambda}$. Moreover, the correlation also includes a phase $e^{ip\cdot
x}$ which becomes $\cos(p\cdot x)$ in the case of the valence quark
distribution. Another important takeaway is that $\bar{\Lambda}$ is a Lorentz-
invariant quantity and should be independent of the external momentum.

However, it must be pointed out that the conclusion in Eq. (38) is based on a
rather crude approximation that $f_{1}$ and $f_{2}$ are constant in $v\cdot
p$. In practice, the transition form factors could have a pole at the mass of
a heavy-light meson created by the current $\bar{q}\gamma^{\mu}h_{v}$ or
$\bar{h}_{v}q$, which is different from $m_{H}$ for the intermediate state
$|H(v)\rangle$. As a result, the binding energy $\bar{\Lambda}$ would also be
different. If we take this into account in Eq. (B.1), then the result will
exhibit a more complicated asymptotic behavior at large distance,

$\displaystyle\lim_{|x|\to\infty}f_{p}(p\cdot x,x^{2})$
$\displaystyle\propto{e^{-\bar{\Lambda}|x|}\over|x|^{d}}\ g[p\cdot
x,\cos(p\cdot x),\sin(p\cdot x)]\,,$ (39)

where the decay constant $\bar{\Lambda}$ should vary among the different
binding energies for the heavy-light mesons, which is similar to the
observation in Ref. [56], and $g$ is a function that can have both oscillating
and non-oscillating dependence on $p\cdot x$. For large enough $|x|$, the
exponential decay should suppress the correlation and make it or its extremes
decrease monotonically in magnitude.

Note that after we match the hybrid scheme matrix elements to
${\overline{\mathrm{MS}}}$, the renormalon ambiguity in the Wilson line mass,
$m_{0}^{{\overline{\mathrm{MS}}}}$, is subtracted out, so the matched result
should exhibit an asymptotic behavior that goes as
$e^{-(\bar{\Lambda}-m_{0}^{{\overline{\mathrm{MS}}}})|z|}$ at large $z$.
Therefore, the sign of $(\bar{\Lambda}-m_{0}^{{\overline{\mathrm{MS}}}})$
becomes crucial in determining whether it is exponentially decaying or
growing.

In QCD sum rule calculations, the result is $\bar{\Lambda}=0.4-0.6$ GeV from
phenomenology, while $m^{\overline{\mathrm{MS}}}_{0}$ is expected to be
$0.1-0.2$ GeV [57], so $\bar{\Lambda}-m^{\overline{\mathrm{MS}}}_{0}=0.2-0.5$
GeV. Since the quarks have heavier-than-physical masses in our lattice
calculation, one should expect a larger $\bar{\Lambda}$, so it is very likely
that $\bar{\Lambda}-m^{\overline{\mathrm{MS}}}_{0}$ still remains positive.
After all, this can be always put to test on the $P^{z}=0$ matrix elements
since $\bar{\Lambda}-m^{\overline{\mathrm{MS}}}_{0}$ is a Lorentz-invariant
quantity.

###  B.2 Extrapolation and FT

If $z_{L}$ is large enough for the correlation $\tilde{h}(z)$ to reach the
asymptotic region, then an extrapolation that encodes the exponential decay
behavior we derived in App. B.1 should lead to reliable FT for moderate-to-
large $x$. To be more precise, there is a rigorous upper bound for the
uncertainty of FT which decreases with $x$.

To prove the above statement, let us consider extrapolation based on the
general model

$\displaystyle\tilde{h}(\lambda)$
$\displaystyle=e^{-c|\lambda-\lambda_{L}|}g(\lambda)\,,$ (40)

where $g(\lambda_{L})=\tilde{h}(\lambda_{L})$, and $c=m_{\rm eff}/P^{z}$ with
$m_{\rm eff}$ being the effective mass for the exponential decay. Motivated by
QCD sum rule results, we expect $m_{\rm eff}\sim 0.2-0.5$ GeV, which can be
larger since we have used heavier-than-physical quark masses. Therefore, for
$P^{z}\sim 2.0$ GeV in the current work, we should have $c\sim 0.10-0.25$ or
higher.

Now let us compare two extrapolations $h_{1}$ and $h_{2}$ with different
$g_{1}$ and $g_{2}$. The difference between the two extrapolations,

$\displaystyle\delta\tilde{h}(\lambda)$
$\displaystyle\equiv\tilde{h}_{1}(\lambda)-\tilde{h}_{2}(\lambda)\,,$ (41)

should satisfy $\delta\tilde{h}(\lambda_{L})=0$ and
$\delta\tilde{h}(\infty)=0$. The difference in the FT with extrapolation is
therefore

$\displaystyle\delta\tilde{f}(x)$
$\displaystyle=\int_{\lambda_{L}}^{\infty}{d\lambda\over\pi}\delta\tilde{h}(\lambda)\cos(x\lambda)\,.$
(42)

If we can approximate $\delta\tilde{h}(\lambda)$ as a flat curve within one
period of the oscillatory function $\cos(x\lambda)$, then the integral in that
region vanishes. This condition can be satisfied if
$|\delta\tilde{h}^{\prime}(\lambda)|\ll x$, which should be reached very
quickly due to the exponential suppression at large $\lambda$. For each $x$,
there should be a minimal integer $N_{x}$ which satisfies
$|\delta\tilde{h}^{\prime}(\lambda_{L}+N_{x}2\pi/x)|\ll x$, so that we can
approximate $\delta\tilde{f}(x)$ as

$\displaystyle\delta\tilde{f}(x)$
$\displaystyle\approx\int_{\lambda_{L}}^{\lambda_{L}+N_{x}{2\pi\over
x}}{d\lambda\over\pi}\delta\tilde{h}(\lambda)\cos(x\lambda)\,.$ (43)

Since $\delta\tilde{h}(\lambda_{L})=0$ and $\delta\tilde{h}(\infty)=0$, there
must be at least one extremum of $\delta\tilde{h}(\lambda)$ for
$\lambda_{L}<\lambda<\infty$, so we have the inequality

$\displaystyle|\delta\tilde{f}(x)|$
$\displaystyle<\int_{\lambda_{L}}^{\lambda_{L}+N_{x}{2\pi\over
x}}{d\lambda\over\pi}\ |\delta\tilde{h}(\lambda)||\cos(x\lambda)|$
$\displaystyle<N_{x}|\delta\tilde{h}(\lambda)|_{\rm
max}\int_{\lambda_{L}}^{\lambda_{L}+{2\pi\over
x}}{d\lambda\over\pi}|\cos(x\lambda)|$
$\displaystyle={4N_{x}|\delta\tilde{h}(\lambda)|_{\rm max}\over\pi
x}\lesssim{4N_{x}|\tilde{h}(\lambda_{L})|\over\pi x}\,.$ (44)

According to our estimate of $\bar{\Lambda}-m_{0}^{{\overline{\mathrm{MS}}}}$,
$c\gtrsim 0.1$ at $P^{z}\sim 2$ GeV, so

$\displaystyle e^{-cN_{x}(2\pi)/x}\lesssim e^{-0.6N_{x}/x}\,,$ (45)

and $N_{x}\sim{\cal O}(1)$ should be sufficient to satisfy
$|\delta\tilde{h}^{\prime}(\lambda_{L}+N_{x}2\pi/x)|\ll x$ with $0<x<1$.
Therefore, in Eq. (B.2) we demonstrate that there is an upper bound for the
model uncertainty in the FT with exponential extrapolation, which decreases in
$x$. The error is also proportional to $|\delta\tilde{h}(\lambda)|_{\rm max}$
which can be much smaller than $|\tilde{h}(\lambda_{L})|$ that is already
close to zero. If $h(\lambda_{L})=0.1$, $|\delta\tilde{h}(\lambda)|_{\rm
max}=0.05$, and $N_{x}=1$, then we have

$\displaystyle|\delta\tilde{f}(x)|<{0.07\over x}\,,$ (46)

which is less than 0.15 at $x=0.5$ and around 15% of the central value of the
qPDF as we obtain below. It is worth pointing out that our estimate of the
upper bound in Eq. (B.2) can be highly overestimated, as
$\delta\tilde{h}(\lambda)$ has an oscillation from $\cos(\lambda)$ and
$\sin(\lambda)$ which are out of pace with $\cos(x\lambda)$ for $0<x<1$, and
$|\delta\tilde{h}(\lambda)|_{\rm max}$ could be much smaller than
$|\tilde{h}(\lambda_{L})|$ and at a sharp peak within
$\lambda_{L}<\lambda<\lambda_{L}+N_{x}2\pi/x$.

Therefore, the FT with exponential extrapolation is under control for moderate
and large $x$. When $\tilde{h}(\lambda_{L})$ is small enough, the model
uncertainty from the extrapolation can be controlled to be much smaller than
the other systematic uncertainties which are about $10\%-20\%$ in this work.

It is worth to compare with the extrapolation error when the correlation
function decreases algebraically as $1/|\lambda|^{d}$, which corresponds to
the generic model

$\displaystyle\tilde{h}(\lambda)$
$\displaystyle=\left({\lambda_{L}\over\lambda}\right)^{d}g(\lambda)\,.$ (47)

Suppose we truncate at $\lambda_{L}=10$, then

$\displaystyle\left({\lambda_{L}\over\lambda_{L}+N_{x}2\pi/x}\right)^{d}$
$\displaystyle\sim(1+0.6N_{x}/x)^{-d}\,.$ (48)

The power $d$ is related to the small-$x$ behavior of the PDF. If we
parameterize the PDF as $\sim x^{a}(1-x)^{b}$, then with LO matching one can
derive that $d={\rm min}\\{1+a,1+b\\}$ [27], which is ${\cal O}(1)$
empirically. Therefore, it will take $N_{x}\gg 1$ for the factor in Eq. (48)
to decrease sufficiently to satisfy the condition
$|\delta\tilde{h}^{\prime}(\lambda_{L}+N_{x}2\pi/x)|\ll x$. As a result, the
uncertainty in the FT is of orders of magnitude larger than that of
extrapolation with exponential decay.

To test our claim of controlled FT error with exponential decay, we choose a
particular model

$\displaystyle\tilde{h}(\lambda)$
$\displaystyle=\tilde{h}(\lambda_{L})\left({\lambda_{L}\over\lambda}\right)^{d}e^{-c|\lambda-\lambda_{L}|}\,.$
(49)

Suppose that the extrapolation is done at $\lambda_{L}=10$ with
$\tilde{h}(\lambda_{L})=0.15$, and the parameters $c$ and $d$ are fitted with
errors $\delta c$ and $\delta d$, then we analytically FT the extrapolated
result to the $x$-space, and calculate its error using

$\displaystyle\delta\tilde{f}(x,c,d)$
$\displaystyle=\sqrt{\left({\partial\tilde{f}\over\partial c}\right)^{2}\delta
c^{2}+\left({\partial\tilde{f}\over\partial d}\right)^{2}\delta d^{2}}\,.$
(50)

Figure 9: Estimate of error in the FT with extroplation using the model in Eq.
(49).

In Fig. 9, we plot the extrapolation error against $x$. We have chosen
different central values of the parameters $c$ and $d$ and fairly large
uncertainties in them. The parameter $d$ cannot have a large negative value,
otherwise it would make $\tilde{h}(\lambda)$ grow beyond $\lambda_{L}$. In
most of the scenarios considered, the error is $\lesssim 0.1$ for $x>0.1$. As
we shall see below, the actual extrapolation error is much smaller than this
estimate and thus negligible when compared to the other systematic errors.

In the following, we perform the extrapolation with four different models. The
extrapolation is carried out on each bootstrap sample by a minimal-square fit.
For each $P^{z}$, we truncate $\tilde{h}(z)$ at the largest $z$, $z_{>0}$,
where the central value of $\tilde{h}(z)$ remains positive, and choose $z_{\rm
max}=\\{z_{>0}-2a,z_{>0}-a,z_{>0}\\}$ to estimate the truncation error. The
range of $z$ used to fit the parameters is $z_{\rm min}\leq z\leq z_{\rm max}$
where $z_{\rm min}$ satisfies $\tilde{h}(z_{\rm min})<0.2$. The continuty
condition between data and model was imposed in the middle point of the fit
range, namely $z_{L}$, which is listed in Table 2. The extrapolation models
are:

| $z_{L}/a$
---|---
$n_{z}$ |  $a=0.04$ fm |  $a=0.06$ fm
$1$ |  $\\{29,30,31\\}$ |  N/A
$2$ |  $\\{26,27,28\\}$ |  $\\{19,20,21\\}$
$3$ |  $\\{19,20,21\\}$ |  $\\{16,17,18\\}$
$4$ |  $\\{24,25,26\\}$ |  $\\{14,15,16\\}$
$5$ |  $\\{21,22,23\\}$ |  $\\{15,16,17\\}$
Table 2: Choices of $z_{L}$ for the extrapolations.

#### Exponential decay model,

or “model-exp”. The model for extrapolation is

$\displaystyle A{e^{-m_{\rm eff}|z|}\over|\lambda|^{d}}\,.$ (51)

We have tried to fit $m_{\rm eff}$ from the same range of $z$ for $P^{z}=0$
matrix elements with a similar form, $Ae^{-m_{\rm eff}|z|}/|z|^{d}$, and found
that $m_{\rm eff}$ is around $0.1$ GeV, about the same scale as the
phenomenological estimate. For the $P^{z}\neq 0$ matrix elements, we do not
fix $m_{\rm eff}$, but constrain it with a prior $m_{\rm eff}\geq m_{\rm
min}$. To test the dependence on this prior condition, we have set $m_{\rm
min}=\\{0,0.1,0.2\\}$ GeV. Besides, we also impose $A>0$ and $d>0$ to ensure
that the extrapolated result is positive and decreases in $\lambda$.

#### Power-law decay model,

or “model-pow”. The model is defined by setting $m_{\rm eff}=0$ in model-exp.
As the $P^{z}\to\infty$ limit of model-exp, model-pow can be used to give a
coarse estimate of the significance of higher-twist effects, although its FT
error is not well under control as we discussed above. We impose the
conditions $A,d>0$ so that the fitted results decrease to zero as
$\lambda\to\infty$.

#### Two-parameter model with exponential decay,

or “model-2p-exp”. As we can see from Fig. 1, the matrix elements at
$\lambda_{L}\sim 6-10$ do not show a clear exponential decay, although they
can be fitted by the latter with $\chi^{2}/d.o.f<1$ due to the large errors.
This may indicate that there is oscillation in $\tilde{h}(\lambda)$. To
incorporate such dependence, we ignore the higher-twist contributions and
assume that the qPDF is parameterized as

$\displaystyle f_{v}(x;a,b)$
$\displaystyle={\Gamma(2+a+b)\over\Gamma(1+a)\Gamma(1+b)}|x|^{a}(1-|x|)^{b}$
$\displaystyle\qquad\times\theta(|x|)\theta(1-|x|)\,.$ (52)

By doing an inverse FT into the $\lambda$-space, the asymptotic form of
$\tilde{h}_{\rm 2p}(\lambda)$ at large $\lambda$ reads,

$\displaystyle\tilde{h}_{\rm 2p}(\lambda)$ $\displaystyle=A\ {\rm
Re}\left[\frac{\Gamma(1+a)}{(-i|\lambda|)^{a+1}}+e^{i\lambda}\frac{\Gamma(1+b)}{(i|\lambda|)^{b+1}}\right]\,.$
(53)

Then we multiply $\tilde{h}_{\rm 2p}(\lambda)$ with an exponential decay
factor as our model for extrapolation,

$\displaystyle\tilde{h}_{\text{2p-exp}}$ $\displaystyle=\tilde{h}_{\rm
2p}(\lambda)e^{-m_{\rm eff}(z-z_{L})}\,.$ (54)

#### Two-parameter model,

or “model-2p”. Again, we ignore the exponential decay and use $\tilde{h}_{\rm
2p}$ as the extrapolation model, which can help us estimate the significance
of higher-twist effects.

Figure 10: FT with different $z_{L}$ for model-exp extrapolation (with prior
$m_{\rm eff}>0.1$ GeV) of the NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$ at $z_{S}=0.24$ fm.

In Fig. 10 we compare the FT with different $z_{L}$ for extrapolation with
model-exp and condition $m_{\rm eff}>0.1$ GeV. Except for very small $x$, the
results are consistent, and those at smaller $z_{L}$ have smaller errors
because the error of the matrix element grows with $z$. Therefore, for the
rest of our analysis, we simply use the largest $z_{L}$ for each $P^{z}$.

Figure 11: Extrapolation with different models for the NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$. At $P^{z}=1.94$ GeV, we have
added the comparison with the 2p-exp and 2p models.

In Fig. 11 we show the extrapolations with different models, which have
noticeable differences at $\lambda>\lambda_{L}$. In Fig. 12 we compare the FT
with different extrapolation models as well as with the discrete FT (DFT). As
we can see, the DFT introduces unphysical oscillation in the qPDF which is due
to the truncation of $\tilde{h}(\lambda)$ at $\lambda_{L}$. In contrast, the
extrapolations are free of such oscillation, and different models yield
consistent qPDFs at moderate and large $x$, thought they differ significantly
at small $x$. We notice that the qPDF from model-2p extrapolation still has
slight oscillations despite its agreement with the others, because the
extrapolated result decays too slowly at $\lambda>\lambda_{L}$. As expected,
the models with exponential decay lead to regular qPDFs at $x=0$, whereas
model-pow and model-2p give divergent qPDFs as $x\to 0$.

Figure 12: Comparison of DFT and FT with different extrapolation models for
the NNLO-matched $\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$ at $z_{S}=0.24$
fm. At $P^{z}=1.94$ GeV, we have added the comparison with the 2p-exp and 2p
models.

Based on the above results, we use model-exp with $m_{\rm eff}>0.1$ GeV for
the FT in our following analysis. To have a coarse estimate of the
uncertainties from extrapolation model and higher-twist contributions, we look
into the difference between final PDFs matched from qPDFs with model-exp and
model-pow extrapolations.

Recall that although the hybrid-scheme matrix elements
$\tilde{h}(\lambda,\lambda_{S},P^{z})$ should be RG invariant, they can still
depend on $\mu$ due to the fixed-order Wilson coefficients used in the
matching between lattice and ${\overline{\mathrm{MS}}}$ schemes. In Fig. 13,
we compare the qPDFs which are FTs of
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$ obtained at $a=0.04$ fm with
$C_{0}^{\rm NLO}$ and $C_{0}^{\rm NNLO}$. We choose $\mu=1.0$ GeV for
$C_{0}^{\rm NLO}$ and $\mu=2.0$ GeV for $C_{0}^{\rm NNLO}$ as the ansatz in
Eq. (15) appear to best describe the lattice matrix elements according to Fig.
8 at these scales. The results are almost identical to each other, which shows
that the renormalon-inspired model with fixed-order Wilson coefficient can
indeed describe the data within a specific window of $\mu$. At NLO, smaller
$\mu$ is favored as $\alpha_{s}(\mu)$ is larger so that the renormalon effects
become important at lower orders. In Fig. 14 we show the $\mu$-dependence of
the qPDFs from NLO- and NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$. As one can see, the results have
mild dependence on $\mu$ which becomes more significant at lower scales.
Therefore, the uncertainty from scale variation will also be larger in this
region.

Figure 13: Comparison of the qPDF with model-exp extrapolation (with $m_{\rm
eff}>0.1$ GeV) of the NLO- and NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$ at $z_{S}=0.24$ fm and
$z_{L}=26a$. The choices of $\mu$ are based on where the renormalon model best
describes the matrix elements.

Figure 14: Comparison of the qPDF at different $\mu$ with model-exp
extrapolation of the NLO- and NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$.

##  Appendix C Perturbative matching

In this section we perform the perturbative matching to the qPDF. Recall that
Eq. (Lattice QCD Determination of the Bjorken-$x$ Dependence of Parton
Distribution Functions at Next-to-next-to-leading Order) relates the qPDF to
the PDF,

$\displaystyle f_{v}(x,\mu)=$
$\displaystyle\int_{-\infty}^{\infty}\frac{dy}{|y|}\
C^{-1}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)\tilde{f}_{v}(y,z_{S},P^{z})$
$\displaystyle\qquad+{\cal
O}\Big{(}\frac{\Lambda_{\text{QCD}}^{2}}{(xP^{z})^{2}},\frac{\Lambda_{\text{QCD}}^{2}}{((1-x)P^{z})^{2}}\Big{)}\,.$
(55)

The matching kernel $C$ can be expanded to $O(\alpha_{s})$ as

$\displaystyle C\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)$
$\displaystyle=\delta\left(\frac{x}{y}-1\right)+\alpha_{s}C^{(1)}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)$
$\displaystyle\quad+\alpha_{s}^{2}C^{(2)}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)+{\cal
O}(\alpha_{s}^{3})\,.$ (56)

The inverse matching kernel $C^{-1}$ can obtained by solving

$\displaystyle\int{dz\over|z|}C^{-1}\\!\left(\frac{x}{z},\frac{\mu}{zP^{z}},|z|\lambda_{S}\right)C\\!\left(\frac{z}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)$
$\displaystyle=\delta\big{(}{x\over y}-1\big{)}$ (57)

order by order in $\alpha_{s}$ [60], and the result is

$\displaystyle
C^{-1}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)$
$\displaystyle=\delta\left(\frac{x}{y}-1\right)-\alpha_{s}C^{(1)}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)$
$\displaystyle\quad+\alpha_{s}^{2}\int{dz\over|z|}C^{(1)}\\!\left(\frac{x}{z},\frac{\mu}{zP^{z}},|z|\lambda_{S}\right)C^{(1)}\\!\left(\frac{z}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)$
$\displaystyle\quad-\alpha_{s}^{2}C^{(2)}\\!\left(\frac{x}{y},\frac{\mu}{yP^{z}},|y|\lambda_{S}\right)+{\cal
O}(\alpha_{s}^{3})\,.$ (58)

It has been shown in Ref. [60] that the inverse matching coefficient satisfies
the correct RG and $P^{z}$-evolution equations.

###  C.1 Numerical implementation of matching

Since in the asymptotic regions,

$\displaystyle\lim_{y\to\infty}C\Big{(}{x\over y}\Big{)}\to{\rm\
finite}\,,\quad\quad\lim_{y\to 0}C\Big{(}{x\over y}\Big{)}\propto{y^{2}\over
x^{2}}\,,$ (59)

and

$\displaystyle C\left({x\over y}\right)$ $\displaystyle\equiv
C_{r}\left({x\over y}\right)-\delta\left({x\over
y}-1\right)\int_{-\infty}^{\infty}dy^{\prime}\ C_{r}(y^{\prime})$ (60)

is a plus function (with “$r$” denotes the $x\neq y$ part) that regulates the
singularity at $y=x$, the convolution integral in Eq. (C) is convergent and
insensitive to the cutoffs for $y\to 0,x,\infty$, as long as the qPDF is
integrable. Therefore, we are able to evaluate the integral numerically within
a finite range of $y$ with a target precision.

Figure 15: Comparison of matrix multiplication to direct numerical integration
for the NLO matching correction to one qPDF sample.

The numerical integration in Eq. (Lattice QCD Determination of the Bjorken-$x$
Dependence of Parton Distribution Functions at Next-to-next-to-leading Order)
is time consuming, especially when we have to perform the matching for the
qPDF on each bootstrap sample. Therefore, to speed up the matching procedure,
we discretize the integral in Eq. (Lattice QCD Determination of the
Bjorken-$x$ Dependence of Parton Distribution Functions at Next-to-next-to-
leading Order) and reexpress it as multiplication of a matching matrix and the
qPDF vector. In our implementation, our integration domain is $-2.0<y<2.0$
discretized with a step size $\delta y=0.001$. Since the qPDF falls very close
to zero at $|y|=2.0$, the corresponding uncertainty is negligible as we have
varied the truncation point. Note that the matching coefficient is a plus
function, the step size $\delta y$ also serves as a soft cutoff for the
singularity near $|x/y|=1$ in the plus functions. To test how well the matrix
multiplication can reproduce the exact numerical intergration, we compare the
NLO corrections to the qPDF from one bootstrap sample using the two methods in
Fig. 15. With our current step size, the results are almost indentical for $x$
as small as $0.01$.

Moreover, to test the reliability of our inverse matching coefficient, which
is obtained through expansion in $\alpha_{s}$, we compare it to direct matrix
inversion. To be specific, we construct a square matching matrix $C$ in $x$
and $y$ with $x,y\in[-2,2]$, which is asymmetric but has dominant diagonal
elements, and then invert it to obtain the inverse matching matrix $C^{-1}$.
At small $\alpha_{s}$, the matrix $C$ can be schematically expressed as

$\displaystyle C$ $\displaystyle={\cal I}+{\cal E}\,,$ (61)

where ${\cal I}$ is an identity matrix, whereas ${\cal E}$ is ${\cal
O}(\alpha_{s})$, so that its inverse can be expanded as

$\displaystyle C^{-1}$ $\displaystyle={\cal I}-{\cal E}+{\cal E}^{2}-{\cal
E}^{3}+\ldots\,.$ (62)

(a)

(b)

Figure 16: (a) Comparison of the NLO matching correction to the qPDF with
matrix inversion and the expansion in Eq. (62) to order $n$. (b) Comparison of
NLO and NNLO matching corrections to the qPDF from direct matrix inversion and
the $\alpha_{s}$-expansion in Eq. (C).

In Fig. 16a we first test the convergence of the solution in Eq. (62) for the
NLO matching matrix. By expanding the solution to order $n$, we calculate the
NLO matching correction to a qPDF sample, and then compare it to the result
from direct matrix inversion. Since our main purpose is to compare the two
inversion methods, we increase the step size to $\delta y=0.01$ to reduce the
computing time regardless the accuracy of numerical integration. We find that
by increasing $n$, the expansion method gradually coverges to direct
inversion, as expected. Of course, in perturbation theory, we should calculate
the matching coefficient to $n$-loop accuracy for consistency, for
$\alpha_{s}$ is the actual power-counting parameter.

In Fig. 16b we compare the NLO and NNLO matching corrections to a qPDF sample
using direct matrix inversion and the $\alpha_{s}$-expansion methods. The
results are basically consistent with each other for almost the entire range
of $x\in(0,1)$, except for small deviations. This is because direct matrix
inversion includes all-order terms in $\alpha_{s}$, and the deviations reflect
the size of higher-order effects, whose smallness shows that the perturbation
series is convergent. With our current two-loop accuracy, we adopt the
$\alpha_{s}$-expansion method.

###  C.2 Perturbative convergence

In Fig. 17 we show the matched results for the PDF from the qPDF obtained from
model-exp extrapolation (with $m_{\rm eff}>0.1$ GeV) of the NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$. As one can see, the NNLO
correction is generally smaller than the NLO correction for moderate $x$,
which indicates good perturbative convergence. Near the end-point regions, the
NLO and NNLO corrections become larger than 50%, which suggests that higher-
order corrections or resummation effects become important.

Figure 17: Upper row: the PDFs from NLO and NNLO matching corrections are
compared to the qPDF (or LO PDF), which is obtained from model-exp (with
$m_{\rm eff}>0.1$ GeV) extrapolation of the NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$. Lower row: the ratio of NLO and
NNLO corrections to the qPDF. Figure 18: Comparison of the PDFs at different
$\mu$ obtained from the NLO- and NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$.

To see whether the NNLO matching reduces the uncertainty from scale variation,
we match qPDFs at different $\mu$ to the corresponding PDFs, and then use
DGLAP equation to evolve the results to $\mu=2.0$ GeV. We use NLO matching
coefficient and LO DGLAP evolution kernel for the qPDF obtained from the NLO-
matched $\tilde{h}(\lambda,\lambda_{s},P^{z},\mu,a)$, and NNLO matching
coefficient and NLO DGLAP evolution kernel for the qPDF obtained from the
NNLO-matched $\tilde{h}(\lambda,\lambda_{s},P^{z},\mu,a)$. The NLO DGLAP
evolution formula takes the following form,

$\displaystyle f_{v}(x,\mu)$ $\displaystyle=f_{v}(x,\mu_{0})$ (63)
$\displaystyle+{\alpha_{s}(\mu_{0})t\over
2\pi}\int_{x}^{1}{dy\over|y|}P_{qq}^{(0)}\left({x\over
y}\right)f_{v}(y,\mu_{0})$ $\displaystyle+\left({\alpha_{s}(\mu_{0})t\over
2\pi}\right)^{2}\int_{x}^{1}{dy\over|y|}\left[P_{qq}^{V(1)}+{1\over
2}P_{qq}^{(0)}\otimes P_{qq}^{(0)}\right.$
$\displaystyle\qquad\qquad\qquad\left.-{\beta_{0}\over
2}P_{qq}^{(0)}\right]\left({x\over y}\right)f_{v}(y,\mu_{0})\,,$ (64)

where $t=\ln(\mu^{2}/\mu_{0}^{2})$, $\beta_{0}=(11C_{A}-2n_{f})/6$,
$P_{qq}^{(0)}$ is the LO splitting kernel, and $P_{qq}^{V(1)}$ is the NLO
splitting kernel [77] for the valence quark PDF.

Figure 19: Comparison of the PDFs obtained from NNLO matching of the qPDFs at
different $\mu$ and NLO DGLAP evolution to $\mu=2.0$ GeV.

Since there are only a few common $\mu$ values for the NLO- and NNLO- matched
$\tilde{h}(\lambda,\lambda_{s},P^{z},\mu,a)$, we choose $\mu=1.4$ and $2.0$
GeV for our comparison. In Fig. 18 we show the scale variation of the PDFs
from NLO and NNLO matching, where only the central values are plotted for our
purpose. As one can see, the NNLO matching correction significantly reduces
the uncertainty for $x\lesssim 0.4$ at NLO, while for $x\gtrsim 0.4$ the NNLO
uncertainty band is still about a factor of one half of the NLO case.
Therefore, the NNLO matching does indeed improve the perturbation theory
uncertainty.

Finally, for the NNLO matching we vary $\mu=2.0$ GeV by a factor of $\sqrt{2}$
and $1/\sqrt{2}$, and then use NLO DGLAP equation to evolve the matched
results to $\mu=2.0$ GeV, whose central vavlues are shown in Fig. 19. As one
can see, there is virtually no difference between choosing $\mu=2.0$ and $2.8$
GeV as the factorization scale, but the lower choice of $\mu=1.4$ GeV does
introduce larger uncertainty mainly because $\alpha_{s}$ becomes too large.
Nevertheless, such uncertainty is still quite small compared to the other
systematics.

###  C.3 Dependence on $P^{z}$, $a$ and extrapolation model

Figure 20: The PDFs from NNLO matching of the qPDFs at different $P^{z}$,
which is obtained from model-exp extrapolation of the NNLO-matched
$\tilde{h}(\lambda,\lambda_{S},P^{z},\mu,a)$.

Figure 21: Comparison of the final results from qPDFs obtained by different
model extrapolations for the FT.

In Fig. 20 we show the $P^{z}$-dependence of the PDF with NNLO matching
correction. We find that despite the considerable differences between the
qPDFs at $P^{z}\leq 1.45$ GeV and those at $P^{z}\geq 1.94$ GeV, the matching
corrections bring the final results closer, which shows the effectiveness of
LaMET. Note that the matching drives the qPDF closer to the smaller $x$
region, so the error bands of the PDFs also shrink after matching as they are
contributed from the larger $x$ region. Moreover, we find that the PDFs start
to converge at $P^{z}\geq 1.29$ GeV, which corresponds to a boost factor of
$\sim 4$. As $P^{z}$ increases, the results becomes smaller as $x\to 1$, which
agrees with our expectation that large momentum suppresses the higher-twist
contributions. It is worth mentioning that both the $P^{z}$-dependence and
matching correction appear to be small for $x$ as low as $0.05$, which hints
that the power correction and resummation effects are less severe than our
naive estimate through power counting.

In Fig. 21 we compare the PDFs matched from the qPDFs with model-exp (with
$m_{\rm eff}>0.1$ GeV) and model-pow extrapolations. For $a=0.04$ fm and
$P^{z}=1.94$ GeV, we also added comparison to the model-2p-exp and model-2p
extrapolations. Despite the differences between the qPDFs at small $x$, the
matched results are almost identical even at the smallest $x$ shown in the
plot. Again, this is the outcome of the PDF receiving contributions from the
qPDF at larger $x$ through matching, which suggests that the extrapolation
error can still be under control for $x$ as small as $\sim 0.01$. Note that
the result from model-2p also shows agreement, but it includes slight
oscillations in the $x$-space, because the extrapolated $\tilde{h}(\lambda)$
decays too slowly in the coordinate space. Therefore, in the region where
other systematic errors are under control, the difference between model-exp
and other extrapolations is negligible, and we will use the model-exp
extrapolation to obtain the final results.

##  Appendix D Final results

The central value of our final result is obtained from the qPDF at $a=0.04$
fm, $z_{S}=0.24$ fm, $z_{L}=0.92$ fm, $\mu=2.0$ GeV and $P^{z}=2.42$ GeV with
exponential extrapolation ($m_{\rm eff}>0.1$ GeV) and NNLO matching. The error
from variation of the factorization scale is obtained by repeating the same
procedure for $\mu=1.4$ and $2.8$ GeV and evolving the matched results to
$\mu=2.0$ GeV with the NLO DGLAP equation, as shown in Fig. 19, where let the
error band cover all the data sets from the three different factorization
scales.

In order to obtain a target precision of 10%, we aim to control the relative
${\cal O}(\alpha_{s}^{3})$ matching correction at $\mu=2.0$ GeV be smaller
than 5%. By assuming that the perturbation series grows geometrically, it
means that the relative NLO correction should be less than
${}^{3}\sqrt{5\%}=37\%$ and the relative NNLO correction less than 14%. By
comparing to Fig. 17, it means that we should exclude the regions $x<0.03$ and
$x>0.88$.

Figure 22: Estimate of the size of power correction $\alpha(x)/P_{z}^{2}$
(upper panel) and its relative size to the qPDF $\tilde{f}_{v}(x,P^{z}=2.42\
{\rm GeV})$ (lower panel).

To estimate the size of the power corrections, we fit the PDFs obtained at
$a=0.04$ fm, $P^{z}=\\{1.45,1.94,2.42\\}$ GeV and $a=0.06$ fm,
$P^{z}=\\{1.72,2.15\\}$ GeV to the ansatz $f_{v}(x)+\alpha(x)/P_{z}^{2}$ for
each fixed $x$, and show the size of the power correction term in Fig. 22. At
$P^{z}=2.42$ GeV, we find that the absolute value of the power correction
diverges at very small $x$, as expected, but its relative size
$\alpha(x)/[P_{z}^{2}f_{v}(x)]$ remains finite because the PDF also diverges.
On the contrary, $\alpha(x)/[P_{z}^{2}f_{v}(x)]$ starts to grow as $x\to 1$.
According to our estimate, $\alpha(x)/[P_{z}^{2}f_{v}(x)]\lesssim 0.1$ for
$0.01<x<0.80$ and $\alpha(x)/[P_{z}^{2}f_{v}(x)]\lesssim 0.05$ for
$0.01<x<0.70$. According to Fig. 21, the qPDF from power-law extrapolation
leads to almost identical PDF after the matching correction for $x$ as small
as 0.01. Our explanation is that the matching correction drives the qPDF to
smaller $x$, so the PDF at a given $x$ receives contributions from the
larger-$x$ region of the qPDF which has less $P^{z}$ dependence. Although
there are logarithms of $\mu/(xP^{z})$ in the matching coefficient which
become large at small $x$, they are always multiplied by the DGLAP splitting
function, which when convoluted with the qPDF always drives the result to
smaller $x$, thus the perturbative correction remains small even at $x=0.03$.

Figure 23: Statistical and scale-variation uncertainty of the PDF obtained
from the qPDF at $a=0.04$ fm and $P^{z}=2.42$ GeV.

In Fig. 23 we show uncertainty of the PDF, $\delta f_{v}(x)/f_{v}(x)$, where
$\delta f_{v}(x)$ includes both statistical and scale-variation errors. The
uncertainty is $\leq 20\%$ for $0.01\leq x\leq 0.93$, as $x=0.01$ is the
smallest $x$ that we show in the plot, and $\leq 10\%$ for $0.08\leq x\leq
0.45$.

Therefore, by combining the estimates of power correction, higher-order
perturbative correction, statistical and scale-variation errors, we determine
the PDF at $0.03\lesssim x\lesssim 0.80$ with $\leq 20\%$ uncertainty and at
$0.08\lesssim x\lesssim 0.45$ with $\leq 10\%$ uncertainty, which is shown in
Fig. 4.

## References

  * Accardi  _et al._ [2016] A. Accardi _et al._ , Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Abdul Khalek  _et al._ [2021] R. Abdul Khalek _et al._ , Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, (2021), arXiv:2103.05419 [physics.ins-det] .
  * Lin  _et al._ [2018] H.-W. Lin _et al._ , Parton distributions and lattice QCD calculations: a community white paper, Prog. Part. Nucl. Phys. 100, 107 (2018), arXiv:1711.07916 [hep-ph] .
  * Constantinou  _et al._ [2021] M. Constantinou _et al._ , Parton distributions and lattice-QCD calculations: Toward 3D structure, Prog. Part. Nucl. Phys. 121, 103908 (2021), arXiv:2006.08636 [hep-ph] .
  * Ji [2013] X. Ji, Parton Physics on a Euclidean Lattice, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Ji [2014] X. Ji, Parton Physics from Large-Momentum Effective Field Theory, Sci. China Phys. Mech. Astron. 57, 1407 (2014), arXiv:1404.6680 [hep-ph] .
  * Ji  _et al._ [2021a] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Large-momentum effective theory, Rev. Mod. Phys. 93, 035005 (2021a), arXiv:2004.03543 [hep-ph] .
  * Liu and Dong [1994] K.-F. Liu and S.-J. Dong, Origin of difference between anti-d and anti-u partons in the nucleon, Phys. Rev. Lett. 72, 1790 (1994), arXiv:hep-ph/9306299 .
  * Detmold and Lin [2006] W. Detmold and C. J. D. Lin, Deep-inelastic scattering and the operator product expansion in lattice QCD, Phys. Rev. D 73, 014501 (2006), arXiv:hep-lat/0507007 .
  * Braun and Müller [2008] V. Braun and D. Müller, Exclusive processes in position space and the pion distribution amplitude, Eur. Phys. J. C 55, 349 (2008), arXiv:0709.1348 [hep-ph] .
  * Chambers  _et al._ [2017] A. J. Chambers, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, A. Schiller, K. Somfleth, R. D. Young, and J. M. Zanotti, Nucleon Structure Functions from Operator Product Expansion on the Lattice, Phys. Rev. Lett. 118, 242001 (2017), arXiv:1703.01153 [hep-lat] .
  * Radyushkin [2017] A. V. Radyushkin, Quasi-parton distribution functions, momentum distributions, and pseudo-parton distribution functions, Phys. Rev. D 96, 034025 (2017), arXiv:1705.01488 [hep-ph] .
  * Ma and Qiu [2018a] Y.-Q. Ma and J.-W. Qiu, Exploring Partonic Structure of Hadrons Using ab initio Lattice QCD Calculations, Phys. Rev. Lett. 120, 022003 (2018a), arXiv:1709.03018 [hep-ph] .
  * Alexandrou  _et al._ [2019] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, A. Scapellato, and F. Steffens, Systematic uncertainties in parton distribution functions from lattice QCD simulations at the physical point, Phys. Rev. D 99, 114504 (2019), arXiv:1902.00587 [hep-lat] .
  * Constantinou and Panagopoulos [2017] M. Constantinou and H. Panagopoulos, Perturbative renormalization of quasi-parton distribution functions, Phys. Rev. D 96, 054506 (2017), arXiv:1705.11193 [hep-lat] .
  * Stewart and Zhao [2018] I. W. Stewart and Y. Zhao, Matching the quasiparton distribution in a momentum subtraction scheme, Phys. Rev. D 97, 054512 (2018), arXiv:1709.04933 [hep-ph] .
  * Alexandrou  _et al._ [2017] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos, and F. Steffens, A complete non-perturbative renormalization prescription for quasi-PDFs, Nucl. Phys. B 923, 394 (2017), arXiv:1706.00265 [hep-lat] .
  * Chen  _et al._ [2018] J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Parton distribution function with nonperturbative renormalization from lattice QCD, Phys. Rev. D 97, 014505 (2018), arXiv:1706.01295 [hep-lat] .
  * Orginos  _et al._ [2017] K. Orginos, A. Radyushkin, J. Karpie, and S. Zafeiropoulos, Lattice QCD exploration of parton pseudo-distribution functions, Phys. Rev. D 96, 094503 (2017), arXiv:1706.05373 [hep-ph] .
  * Braun  _et al._ [2019] V. M. Braun, A. Vladimirov, and J.-H. Zhang, Power corrections and renormalons in parton quasidistributions, Phys. Rev. D 99, 014013 (2019), arXiv:1810.00048 [hep-ph] .
  * Li  _et al._ [2021] Z.-Y. Li, Y.-Q. Ma, and J.-W. Qiu, Extraction of Next-to-Next-to-Leading-Order Parton Distribution Functions from Lattice QCD Calculations, Phys. Rev. Lett. 126, 072001 (2021), arXiv:2006.12370 [hep-ph] .
  * Fan  _et al._ [2020] Z. Fan, X. Gao, R. Li, H.-W. Lin, N. Karthik, S. Mukherjee, P. Petreczky, S. Syritsyn, Y.-B. Yang, and R. Zhang, Isovector parton distribution functions of the proton on a superfine lattice, Phys. Rev. D 102, 074504 (2020), arXiv:2005.12015 [hep-lat] .
  * Ji  _et al._ [2017] X. Ji, J.-H. Zhang, and Y. Zhao, More On Large-Momentum Effective Theory Approach to Parton Physics, Nucl. Phys. B 924, 366 (2017), arXiv:1706.07416 [hep-ph] .
  * Radyushkin [2018] A. V. Radyushkin, Quark pseudodistributions at short distances, Phys. Lett. B 781, 433 (2018), arXiv:1710.08813 [hep-ph] .
  * Izubuchi  _et al._ [2018] T. Izubuchi, X. Ji, L. Jin, I. W. Stewart, and Y. Zhao, Factorization Theorem Relating Euclidean and Light-Cone Parton Distributions, Phys. Rev. D 98, 056004 (2018), arXiv:1801.03917 [hep-ph] .
  * Gao  _et al._ [2020] X. Gao, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, S. Syritsyn, and Y. Zhao, Valence parton distribution of the pion from lattice QCD: Approaching the continuum limit, Phys. Rev. D 102, 094513 (2020), arXiv:2007.06590 [hep-lat] .
  * Ji  _et al._ [2021b] X. Ji, Y. Liu, A. Schäfer, W. Wang, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, A Hybrid Renormalization Scheme for Quasi Light-Front Correlations in Large-Momentum Effective Theory, Nucl. Phys. B 964, 115311 (2021b), arXiv:2008.03886 [hep-ph] .
  * Hua  _et al._ [2021] J. Hua, M.-H. Chu, P. Sun, W. Wang, J. Xu, Y.-B. Yang, J.-H. Zhang, and Q.-A. Zhang (Lattice Parton), Distribution Amplitudes of K* and $\phi$ at the Physical Pion Mass from Lattice QCD, Phys. Rev. Lett. 127, 062002 (2021), arXiv:2011.09788 [hep-lat] .
  * Chen  _et al._ [2021] L.-B. Chen, W. Wang, and R. Zhu, Next-to-Next-to-Leading Order Calculation of Quasiparton Distribution Functions, Phys. Rev. Lett. 126, 072002 (2021), arXiv:2006.14825 [hep-ph] .
  * Gluck  _et al._ [1992] M. Gluck, E. Reya, and A. Vogt, Pionic parton distributions, Z. Phys. C 53, 651 (1992).
  * Novikov  _et al._ [2020] I. Novikov _et al._ , Parton Distribution Functions of the Charged Pion Within The xFitter Framework, Phys. Rev. D 102, 014040 (2020), arXiv:2002.02902 [hep-ph] .
  * Barry  _et al._ [2021] P. C. Barry, C.-R. Ji, N. Sato, and W. Melnitchouk (Jefferson Lab Angular Momentum (JAM)), Global QCD Analysis of Pion Parton Distributions with Threshold Resummation, Phys. Rev. Lett. 127, 232001 (2021), arXiv:2108.05822 [hep-ph] .
  * Aicher  _et al._ [2010] M. Aicher, A. Schafer, and W. Vogelsang, Soft-gluon resummation and the valence parton distribution function of the pion, Phys. Rev. Lett. 105, 252003 (2010), arXiv:1009.2481 [hep-ph] .
  * Zhang  _et al._ [2019] J.-H. Zhang, J.-W. Chen, L. Jin, H.-W. Lin, A. Schäfer, and Y. Zhao, First direct lattice-QCD calculation of the $x$-dependence of the pion parton distribution function, Phys. Rev. D 100, 034505 (2019), arXiv:1804.01483 [hep-lat] .
  * Sufian  _et al._ [2019] R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu, and D. G. Richards, Pion Valence Quark Distribution from Matrix Element Calculated in Lattice QCD, Phys. Rev. D 99, 074507 (2019), arXiv:1901.03921 [hep-lat] .
  * Izubuchi  _et al._ [2019] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, and S. Syritsyn, Valence parton distribution function of pion from fine lattice, Phys. Rev. D 100, 034516 (2019), arXiv:1905.06349 [hep-lat] .
  * Joó  _et al._ [2019] B. Joó, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, R. S. Sufian, and S. Zafeiropoulos, Pion valence structure from Ioffe-time parton pseudodistribution functions, Phys. Rev. D 100, 114512 (2019), arXiv:1909.08517 [hep-lat] .
  * Sufian  _et al._ [2020] R. S. Sufian, C. Egerer, J. Karpie, R. G. Edwards, B. Joó, Y.-Q. Ma, K. Orginos, J.-W. Qiu, and D. G. Richards, Pion Valence Quark Distribution from Current-Current Correlation in Lattice QCD, Phys. Rev. D 102, 054508 (2020), arXiv:2001.04960 [hep-lat] .
  * Lin  _et al._ [2021] H.-W. Lin, J.-W. Chen, Z. Fan, J.-H. Zhang, and R. Zhang, Valence-Quark Distribution of the Kaon and Pion from Lattice QCD, Phys. Rev. D 103, 014516 (2021), arXiv:2003.14128 [hep-lat] .
  * Gao  _et al._ [2021] X. Gao, K. Lee, S. Mukherjee, C. Shugert, and Y. Zhao, Origin and resummation of threshold logarithms in the lattice QCD calculations of PDFs, Phys. Rev. D 103, 094504 (2021), arXiv:2102.01101 [hep-ph] .
  * Bazavov  _et al._ [2014] A. Bazavov _et al._ (HotQCD), Equation of state in ( 2+1 )-flavor QCD, Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Follana  _et al._ [2007] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), Highly improved staggered quarks on the lattice, with applications to charm physics, Phys. Rev. D 75, 054502 (2007), arXiv:hep-lat/0610092 .
  * Hasenfratz and Knechtli [2001] A. Hasenfratz and F. Knechtli, Flavor symmetry and the static potential with hypercubic blocking, Phys. Rev. D 64, 034504 (2001), arXiv:hep-lat/0103029 .
  * Ji  _et al._ [2018] X. Ji, J.-H. Zhang, and Y. Zhao, Renormalization in Large Momentum Effective Theory of Parton Physics, Phys. Rev. Lett. 120, 112001 (2018), arXiv:1706.08962 [hep-ph] .
  * Ishikawa  _et al._ [2017] T. Ishikawa, Y.-Q. Ma, J.-W. Qiu, and S. Yoshida, Renormalizability of quasiparton distribution functions, Phys. Rev. D 96, 094019 (2017), arXiv:1707.03107 [hep-ph] .
  * Green  _et al._ [2018] J. Green, K. Jansen, and F. Steffens, Nonperturbative Renormalization of Nonlocal Quark Bilinears for Parton Quasidistribution Functions on the Lattice Using an Auxiliary Field, Phys. Rev. Lett. 121, 022004 (2018), arXiv:1707.07152 [hep-lat] .
  * Bauer  _et al._ [2012] C. Bauer, G. S. Bali, and A. Pineda, Compelling Evidence of Renormalons in QCD from High Order Perturbative Expansions, Phys. Rev. Lett. 108, 242002 (2012), arXiv:1111.3946 [hep-ph] .
  * Zhang  _et al._ [2017] J.-H. Zhang, J.-W. Chen, X. Ji, L. Jin, and H.-W. Lin, Pion Distribution Amplitude from Lattice QCD, Phys. Rev. D 95, 094514 (2017), arXiv:1702.00008 [hep-lat] .
  * Green  _et al._ [2020] J. R. Green, K. Jansen, and F. Steffens, Improvement, generalization, and scheme conversion of Wilson-line operators on the lattice in the auxiliary field approach, Phys. Rev. D 101, 074509 (2020), arXiv:2002.09408 [hep-lat] .
  * Alexandrou  _et al._ [2021] C. Alexandrou, K. Cichy, M. Constantinou, J. R. Green, K. Hadjiyiannakou, K. Jansen, F. Manigrasso, A. Scapellato, and F. Steffens, Lattice continuum-limit study of nucleon quasi-PDFs, Phys. Rev. D 103, 094512 (2021), arXiv:2011.00964 [hep-lat] .
  * Huo  _et al._ [2021] Y.-K. Huo _et al._ (Lattice Parton Collaboration (LPC)), Self-renormalization of quasi-light-front correlators on the lattice, Nucl. Phys. B 969, 115443 (2021), arXiv:2103.02965 [hep-lat] .
  * Bazavov  _et al._ [2018a] A. Bazavov, P. Petreczky, and J. H. Weber, Equation of State in 2+1 Flavor QCD at High Temperatures, Phys. Rev. D 97, 014510 (2018a), arXiv:1710.05024 [hep-lat] .
  * Bazavov  _et al._ [2016] A. Bazavov, N. Brambilla, H. T. Ding, P. Petreczky, H. P. Schadler, A. Vairo, and J. H. Weber, Polyakov loop in 2+1 flavor QCD from low to high temperatures, Phys. Rev. D 93, 114502 (2016), arXiv:1603.06637 [hep-lat] .
  * Bazavov  _et al._ [2018b] A. Bazavov, N. Brambilla, P. Petreczky, A. Vairo, and J. H. Weber (TUMQCD), Color screening in (2+1)-flavor QCD, Phys. Rev. D 98, 054511 (2018b), arXiv:1804.10600 [hep-lat] .
  * Petreczky  _et al._ [2021] P. Petreczky, S. Steinbeisser, and J. H. Weber, Chromo-electric screening length in 2+1 flavor QCD, in  _38th International Symposium on Lattice Field Theory_ (2021) arXiv:2112.00788 [hep-lat] .
  * Burkardt  _et al._ [1995] M. Burkardt, J. M. Grandy, and J. W. Negele, Calculation and interpretation of hadron correlation functions in lattice QCD, Annals Phys. 238, 441 (1995), arXiv:hep-lat/9406009 .
  * Beneke and Braun [1994] M. Beneke and V. M. Braun, Heavy quark effective theory beyond perturbation theory: Renormalons, the pole mass and the residual mass term, Nucl. Phys. B 426, 301 (1994), arXiv:hep-ph/9402364 .
  * Xiong  _et al._ [2014] X. Xiong, X. Ji, J.-H. Zhang, and Y. Zhao, One-loop matching for parton distributions: Nonsinglet case, Phys. Rev. D 90, 014051 (2014), arXiv:1310.7471 [hep-ph] .
  * Ma and Qiu [2018b] Y.-Q. Ma and J.-W. Qiu, Extracting Parton Distribution Functions from Lattice QCD Calculations, Phys. Rev. D 98, 074021 (2018b), arXiv:1404.6860 [hep-ph] .
  * [60] Y. Zhao _et al._ , in preparation.
  * Gao [2021] X. Gao,  _Pion Valence Structure and Form Factors from Lattice QCD at Physical Point,_ (Talk given at “9th Workshop of the APS Topical Group on Hadronic Physics”, Virtual, US, 2021).
  * Lin and Zhang [2019] H.-W. Lin and R. Zhang, Lattice finite-volume dependence of the nucleon parton distributions, Phys. Rev. D 100, 074502 (2019).
  * Liu and Chen [2021] W.-Y. Liu and J.-W. Chen, Chiral perturbation for large momentum effective field theory, Phys. Rev. D 104, 054508 (2021), arXiv:2011.13536 [hep-lat] .
  * Braun  _et al._ [2020] V. M. Braun, K. G. Chetyrkin, and B. A. Kniehl, Renormalization of parton quasi-distributions beyond the leading order: spacelike vs. timelike, JHEP 07, 161, arXiv:2004.01043 [hep-ph] .
  * Bali  _et al._ [2013] G. S. Bali, C. Bauer, A. Pineda, and C. Torrero, Perturbative expansion of the energy of static sources at large orders in four-dimensional SU(3) gauge theory, Phys. Rev. D 87, 094517 (2013), arXiv:1303.3279 [hep-lat] .
  * Bali  _et al._ [2014a] G. S. Bali, C. Bauer, and A. Pineda, Perturbative expansion of the plaquette to ${\cal O}(\alpha^{35})$ in four-dimensional SU(3) gauge theory, Phys. Rev. D 89, 054505 (2014a), arXiv:1401.7999 [hep-ph] .
  * Bali  _et al._ [2014b] G. S. Bali, C. Bauer, and A. Pineda, Model-independent determination of the gluon condensate in four-dimensional SU(3) gauge theory, Phys. Rev. Lett. 113, 092001 (2014b), arXiv:1403.6477 [hep-ph] .
  * Bali  _et al._ [2014c] G. S. Bali, S. Collins, B. Gläßle, M. Göckeler, J. Najjar, R. H. Rödl, A. Schäfer, R. W. Schiel, A. Sternbeck, and W. Söldner, The moment $\langle x\rangle_{u-d}$ of the nucleon from $N_{f}=2$ lattice QCD down to nearly physical quark masses, Phys. Rev. D 90, 074510 (2014c), arXiv:1408.6850 [hep-lat] .
  * Pineda [2003] A. Pineda, The Static potential: Lattice versus perturbation theory in a renormalon based approach, J. Phys. G 29, 371 (2003), arXiv:hep-ph/0208031 .
  * Shifman  _et al._ [1979a] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, QCD and Resonance Physics. Theoretical Foundations, Nucl. Phys. B 147, 385 (1979a).
  * Shifman  _et al._ [1979b] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, QCD and Resonance Physics: Applications, Nucl. Phys. B 147, 448 (1979b).
  * Novikov  _et al._ [1984] V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Two-Dimensional Sigma Models: Modeling Nonperturbative Effects of Quantum Chromodynamics, Phys. Rept. 116, 103 (1984).
  * Novikov  _et al._ [1985] V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Wilson’s Operator Expansion: Can It Fail?, Nucl. Phys. B 249, 445 (1985).
  * David [1986] F. David, The Operator Product Expansion and Renormalons: A Comment, Nucl. Phys. B 263, 637 (1986).
  * Petreczky and Weber [2020] P. Petreczky and J. H. Weber, Strong coupling constant from moments of quarkonium correlators revisited, (2020), arXiv:2012.06193 [hep-lat] .
  * Falk  _et al._ [1990] A. F. Falk, H. Georgi, B. Grinstein, and M. B. Wise, Heavy Meson Form-factors From QCD, Nucl. Phys. B 343, 1 (1990).
  * Curci  _et al._ [1980] G. Curci, W. Furmanski, and R. Petronzio, Evolution of Parton Densities Beyond Leading Order: The Nonsinglet Case, Nucl. Phys. B 175, 27 (1980).