Towards High-Precision Parton Distributions From Lattice QCD via Distillation

We apply the Distillation spatial smearing program to the extraction of the unpolarized isovector valence PDF of the nucleon. The improved volume sampling and control of excited-states afforded by distillation leads to a dramatically improved determination of the requisite Ioffe-time Pseudo-distribution (pITD). The impact of higher-twist effects is subsequently explored by extending the Wilson line length present in our non-local operators to one half the spatial extent of the lattice ensemble considered. The valence PDF is extracted by analyzing both the matched Ioffe-time Distribution (ITD), as well as a direct matching of the pITD to the PDF. Through development of a novel prescription to obtain the PDF from the pITD, we establish a concerning deviation of the pITD from the expected DGLAP evolution of the pseudo-PDF. The presence of DGLAP evolution is observed once more following introduction of a discretization term into the PDF extractions. Observance and correction of this discrepancy further highlights the utility of distillation in such structure studies.


Introduction
Elucidation of the dynamical properties of quarks and gluons and their collective emergent phenomena is the principal charge of contemporary hadronic physics and is a central component of the precision phenomenology program at the Large Hadron Collider (LHC). Crucial for the interpretation and prediction of inclusive and semi-inclusive scattering processes, such as deep inelastic scattering (DIS) and semi-inclusive DIS, are the parton distribution functions (PDFs) which emerge in the QCD factorization [1] of such inclusive cross sections. PDFs capture the collinear momentum distribution of partons within a fast moving hadron with p µ = p + , m 2 /2p + , 0 ⊥ , and offer a probabilistic interpretation featured at leadingtwist. Despite complementing the growing hadron tomography efforts both theoretically and at upcoming facilities such as the Electron Ion Collider, determination of PDFs remains a high priority as they are often large sources of hadronic error in collider experiments and thereby affect the precision measurements of an array of Standard Model parameters [2].
The numerical tool of lattice field theory enables the quantitative study of stronglycoupled theories, such as Quantum Chromodynamics (QCD), from first-principles. As PDFs accumulate information on the infrared structure of a hadron, they would seem ideal objects JHEP11(2021)148 with ν ≡ p · z and z 2 the Ioffe-time [34,35] and invariant interval, respectively, of the process. For a fast-moving hadron, the usual unpolarized PDFs are defined with light-cone coordinates where α = +, p α = p + , m 2 h 2p + , 0 ⊥ and z α = (0, z − , 0 ⊥ ). In this scenario M + (p, z) only receives contributions from M p + z − , 0 . Provided the logarithmic singularity that arises for z 2 = 0 is regularized (typically in MS), M p + z − , 0 defines the Ioffe-time distribution (ITD) [35]: and obviates the Fourier transform of (2.1) to momentum space which defines the conventional PDFs. Lorentz invariance implies the ν-dependence of M p + z − , 0 µ 2 can be computed in any frame, and with any choice of {z, α} that may be convenient. A particular choice amenable to calculation with lattice QCD is α = 0, p α = (E, 0 ⊥ , p z ) and z α = (0, 0 ⊥ , z 3 ), which excludes the contamination from the pure higher-twist term N ν, z 2 . The remaining term M ν, z 2 = 0 is deemed the Ioffe-time Pseudo-distribution [29] or pseudo-ITD. In addition to the twist-2 contributions, the pseudo-ITD also contains higher-twist contributions O(z 2 Λ 2 QCD ) that vanish only in the light-cone limit. Furthermore, for all relevant Feynman diagrams [36], the Fourier transform of the pseudo-ITD with respect to ν has support only on the canonical parton momentum fraction interval x ∈ [−1, 1]. The challenge numerically is to extract the leading-twist dependence amongst contributing O z 2 higher-twist terms.
A considerable challenge of using this non-local parton bilinear is the appearance of additional ultraviolet (UV) divergences for space-like separations. Prior to taking the continuum limit, the UV divergences must be regularized and removed. In perturbation theory, such divergences appear first, as power-like z/a [37] terms in the gauge-link selfenergy corrections, with a being a UV regulator. These exponentiate to all orders [38,39] producing the factor Z link (z 3 , a) e −A|z 3 |/a . Second, there are logarithmic ln(−z 2 /a 2 ) UV corrections present both in the link self-energy and in the vertex link corrections [40]. As these unwanted UV divergences are multiplicative [40][41][42][43], and independent of the Ioffe-time, combining into an overall factor Z UV (z 3 , a), we construct the so-called reduced pseudo-ITD [29] M ν, z 2 = M ν, z 2 M (0, z 2 ) . (2.4) Since the rest-frame pseudo-ITD M 0, z 2 has the same UV divergences associated with the gauge link, by forming the reduced pseudo-ITD we cancel the divergent factors Z link (z 3 , a) and Z UV (z 3 , a), thereby ensuring a finite continuum limit. The choice to construct the RG invariant reduced pseudo-ITD with M 0, z 2 is especially straightforward, as M 0, z 2 is simply the bare vector charge Z −1 V in the light-cone limit thereby leaving the OPE unaltered. Motivation for the reduced pseudo-ITD also extends to the mitigation of higher-twist O(z 2 Λ 2 QCD ) effects [44]. An alternative has recently been proposed that makes use of a vacuum matrix element of the space-like parton bilinear [45] (see also ref. [46]).
The factorization is valid in so far as the polynomial corrections B k (ν) z 2 k can be mitigated. We follow in this work the plus-prescription defined by 1

Numerical methods
The pseudo-distribution formalism has been leveraged in several lattice calculations of partonic structure of hadrons, including the valence quark content of the pion [54], and the unpolarized valence quark [44,[55][56][57][58] and recently gluon [59] contents of the nucleon. Even though each calculation makes use of standard spatial and momentum smearing techniques, considerable statistical fluctuations are met for Ioffe-times in excess of ν 5 (and at even smaller values for gluonic matrix elements). Motivated by the demonstrable success of the union of the distillation paradigm and the momentum smearing idea [60], we apply for the first time the distillation spatial smearing program to the extraction of PDFs from lattice QCD. We employ a 349 configuration isotropic clover ensemble generated by the JLab/W&M/LANL collaboration [61], featuring 2 ⊕ 1 quark flavors within a 32 3 × 64 lattice volume. The inverse coupling was set to β = 6.3, from which the lattice spacing a 0.094 fm was obtained from the Wilson-flow scale w 0 [62], which yields a pion mass m π = 358 MeV. The reader is referred to [63,64] for further details of this ensemble, denoted a094m358 and summarized in table 1. A single application of the stout smearing kernel [65] yields a tadpole-improved tree-level clover coefficient c sw near the value determined non-perturbatively from the Schrödinger functional method [61] a posteriori.

Distillation
The Jacobi smearing kernel J σ,nσ (T ) = 1 + σ∇ 2 (T ) nσ [66] featured in a variety of lattice calculations is but one method with which point-like interpolating fields can be spatially and ordering solutions according to the eigenvalue magnitude λ (k) (T ). Forming the outer product of equal-time eigenvectors defines the distillation [67] smearing kernel where R D is the desired distillation space rank and color indices {a, b} are made explicit. Two-point correlation functions formed by Wick-contracting quark fields smeared via (3.2) can be recast into a trace over distinct reusable objects constructed within the distillation space, the so-called elementals and perambulators. The elementals encode the interpolator construction, which in the case of baryons read where D i are covariant derivatives, and S αβγ encode the patterns of subduction of a continuum interpolator across irreducible representations (irreps) of a hypercubic lattice and its associated little groups. Quark propagation between distillation spaces is captured by the perambulators with M the Dirac operator. In the case of three-point correlation functions, an additional computational unit appears. The generalized perambulator, or genprop for short, carries the same external indices as a standard perambulator, but includes an intermediate operator insertion. In the case of the Wilson-line operator specific to this work, the genprop reads where the unpolarized PDFs are selected with the Dirac matrix γ 4 , and the sum over the initial spatial coordinates of the Wilson line y ensures a zero 3-momentum projection. We JHEP11(2021)148 The viability of the pseudo-distribution formalism hinges on the space-like quark bilinear remaining in a perturbative, or short-distance, regime. To then minimize the impact of polynomial-z 2 corrections when mapping the ν-dependence of the reduced pseudo-ITD, the matrix elements (2.1) must be isolated in frames of varying external nucleon boosts. Increasing the 3-momentum of a hadronic interpolator not only leads to exponentially worsening statistical fluctuations, but also reduces the efficacy of spatially-smeared interpolators to overlap onto the desired hadronic states. An improvement program capable of increasing interpolator-state overlaps in boosted frames, now known as momentum smearing was first established in [68]. Momentum smearing effectively shifts the operator-groundstate overlap peak in momentum space; crucially, this shift also improves the overlap of a boosted interpolator onto neighboring excited-states. In anticipation of an increasingly dense nucleon spectrum in boosted frames, we adopt in this work a modification of the distillation paradigm [60] which simultaneously leverages the momentum smearing heuristic with improved control over excited-states. In the interest of self-containment, we now summarize this procedure.
The momentum smearing idea is incorporated within distillation by applying spatiallyvarying phases to a pre-computed eigenvector basis where we designate the modified eigenvectors as phased. As the eigenvector basis already reflects the periodicity of the spatial lattice, the phase factors introduced in this manner are restricted to allowed lattice momenta. This requirement was not found to be limiting, but rather offered broad improvement of momentum space overlaps for a range of nucleon momenta [60]. All distillation components must then be reconstructed on each new modified basis. In the interest of compute cycles and storage, we adopt three eigenvector bases: the pre-computed basis (3.1) for the nucleon at rest and with smallẑ-momenta (|ap z | ≤ 3 [2π/L]), while two additional bases are formed according to

Matrix element isolation
Interpolator construction. The breaking of rotational symmetry by a hypercubic lattice entails baryons appear in lattice calculations according to definite patterns of subduction across the finite number of irreps Λ of the octahedral group double-cover O D h . We elect to use a single spatially local, non-relativistic nucleon interpolating operator constructed according to [69,70]. The paradigm established in [71] is adopted herein to ensure our interpolator transforms irreducibly under the appropriate little group irreps. A forthcoming calculation will leverage an extended basis of interpolators to extract these same matrix elements.
The space-like matrix element (2.1) is isolated in a flavor isovector combination according to the kinematics highlighted in section 2, and requires computation of standard two-point (3.8) and three-point correlation functions featuring the unrenormalized Wilson line operator O where energy eigenstate contributions are denoted, and the nucleon interpolating fields N that are smeared with distillation are separated by a Euclidean time T . The explicit momentum projection built into the genprop (3.5) leads to an overall spatial volume factor V 3 in the forward case. The Wilson line operator is inserted for 0 < τ < T . These correlators and the factorization manifest through distillation are illustrated in figure 1. The spectral representations of (3.8) and (3.9) indicate the desired ground-state matrix element follows from ratios of three-point to two-point correlation functions, which plateau asymptotically for 0 τ T . The contamination from excited-states is reduced further in this calculation by extracting the matrix elements using the Summation method [72,73], whereby the time slice τ on whichO Fit T /a ∈ [4,14] Fit T /a ∈ [6,14] Fit T /a ∈ [8,14] R (p z , z 3 ; T ) (a) Fit T /a ∈ [4,14] Fit T /a ∈ [6,14] Fit T /a ∈ [8,14]  to plateau and multi-state fits. The geometric series resulting from (3.10) depends linearly on the targeted matrix element M 4 (p z , z 3 ), for which we implement the fitting functional We note in practice, the excited-state term O e −∆ET is found to have no impact on our summation fits and is hence omitted from our results.
The two-and three-point functions are computed on four temporal source origins per configuration with T /a ∈ {4, 6, 8, 10, 12, 14} ∼ 0.38-1.32 fm. This number of source-sink separations is chosen to filter out any excited-states that are not captured by the combined effect of distillation and the summation method, and to ensure our linear fits (3.11) do not over fit our data as signal-to-noise problems become unavoidable. We consider nucleon momenta up to |p z | = 6 × [2π/aL] ∼ 2.47 GeV and Wilson line lengths up to z 3 /a = 16, although only z 3 /a ≤ 12 will be subsequently used in our analysis. As stated in table 1, the number of temporal origins and configurations is held fixed for all computed matrix elements regardless of T /a, p z and z 3 . A representative set of R (p 3 , z 3 ; T ) and applied linear fits are shown in figures 2 and 3. Repeating the matrix element extraction for all momenta and displacements, we gather the real/imaginary components of the unpolarized reduced pseudo-ITD in figures 4a and 4b.
Our adoption of distillation in this calculation is justified numerically by making a direct comparison in appendix A with an earlier calculation of the nucleon's unpolarized PDF [58], which instead used spatial Jacobi smearing on the same ensemble. A cost comparison between the two approaches is also presented, making a strong case for distillation in future calculations of hadronic structure. Fit T /a ∈ [4,14] Fit T /a ∈ [6,14] Fit T /a ∈ [8,14] R (p z , z 3 ; T ) (a) Fit T /a ∈ [4,14] Fit T /a ∈ [6,14] Fit T /a ∈ [8,14]

Extraction of unpolarized nucleon PDFs
As PDFs are determined phenomenologically at a factorization scale µ 2 in MS to renormalize the associated collinear divergences, the nucleon unpolarized reduced pseudo-ITD shown in figure 4a and 4b must be matched to a common scale in MS prior to any meaningful comparisons. At one-loop and without loss of generality, negating the sign of the O (α s ) correction and interchanging the ITD and reduced pseudo-ITD in (2.5) one obtains the JHEP11(2021)148 factorization relationship that matches the reduced pseudo-ITD to the ITD: This relationship describes the evolution of each distinct set of M ν, z 2 data at a given z 2 to a common scale µ 2 in MS. Regardless of whether the evolution and matching steps are done separately or in one step, a smooth and continuous description of the reduced pseudo-ITD for each z 2 in the interval [0, ν] is required. It is common in the literature to find polynomials in Ioffe-time fit to each set of distinct z 2 data in order to build M uν, z 2 [56][57][58]74]. Interpolations are also common, and when used have been found to be consistent with polynomial fits [57,58]. A polynomial in ν is perhaps a dubious choice, as it cannot capture the correct limiting behavior of the ITD at large-ν. To understand this, consider a simple nucleon valence PDF ansatz The cosine transform of this ansatz with respect to Ioffe-time is given by with 2 F 3 a generalized hypergeometric function and α, β > −1. In the large Ioffe-time regime (4.3) behaves as For the real component of the ITD to correspond to a valence PDF with a finite sum rule, the ITD must then vanish for asymptotically large-ν (i.e. α > −1). This suggests the usefulness of a smooth polynomial in ν extends only so far as interpolating the discrete pseudo-ITD data, and should not be used as a measure of the moments of the pseudo-PDFs given their divergent behavior at large-ν. This motivates methods to directly extract the PDFs from the reduced pseudo-ITD, thereby obviating the need for a continuous description of M ν, z 2 in order to perform the evolution/matching steps. This will be developed in sections 4.2 and 4.3.
To get a handle on the scale dependence of our data and ground the ensuing discussion, we nonetheless start with a provisional sixth order polynomial fit in Ioffe-time to the reduced pseudo-ITD for constant z 2 : (4.5) The even (odd) powers of the polynomial are applied to jackknife samples of the real (imaginary) component of M ν, z 2 given in figures 4a and 4b. Higher order polynomials were considered, but were found to be unconstrained by the data. With the polynomial fits JHEP11(2021)148 The NLO prefactor α s C F /2π is included in these data, but omitted from the labels for clarity. The convolutions were performed up to z/a = 16, but data for z/a > 9 are generally noisy and not shown.
in hand, we perform the evolution and scheme conversion convolutions (4.1) in a single step. The matched MS scale µ = 2 GeV was chosen, and the strong coupling α s (2 GeV) 0.303 was adopted from LHAPDF6 [75]. The scale µ = 2 GeV corresponds to the reduced pseudo-ITD being evolved to the common scale GeV. On this ensemble a094m358, this common scale then equates to z 2 0 /a 2 0.511. The computed evolution and scheme matching convolutions are depicted in figure 5 and figure 6 for the real and imaginary components, respectively. It is curious the evolution and matching convolutions appear to be nearly equal in magnitude but opposite in sign. This feature of the pseudo-distributions has been observed in independent calculations [56,57] and hints an NNLO matching relation may not be needed. Nonetheless, a future study will explore to what effect the matching relation (2.5), truncated at NLO in this work, can be improved at NNLO [46].
When the scale and scheme conversion are incorporated, we observe in figure 7a and figure 7b a dramatic collapse of the reduced pseudo-ITD onto a common curve for z/a 10. The lack of residual z 2 -dependence is particularly striking in the real component of the ITD, but less so in the imaginary component. This confirms the formation of the reduced ratio (2.4) indeed cancels much of the z 2 -dependence in the pseudo-ITD, with any remaining at small-z 2 ideally described by the coordinate-space DGLAP evolution [29]. The NLO prefactor α s C F /2π is included in these data, but omitted from the labels for clarity. The convolutions were performed up to z/a = 16, but data for z/a > 9 are generally noisy and not shown.

An ill-posed inverse
The Fourier transform relating the x-dependence of the PDF to the ν-dependence of the ITD is an ill-posed inverse problem, as the ITD is computed in a discrete and limited range of Ioffe-time. Regularization at this stage is however numerically cheap and more stable relative to a direct matching of the {x, µ 2 } dependencies of the PDF to the {ν, z 2 } dependencies of the reduced pseudo-ITD. A direct inversion of (4.6) is satisfied by an infinite number of solutions, each of which having little predictive or postdictive credibility. The futility of direct inversions has been demonstrated in a few LCS calculations [57,58], wherein each inversion yielded unstable PDFs with spurious oscillations. The limited range of Ioffe-time accessible to present PDF calculations only compounds the need for refined extraction methods. This inversion problem is shared with other lattice formalisms that rely on QCD factorization, and indeed the global analysis of inclusive/semi-inclusive processes. Although different in character, this problem pervades the quantitative sciences and even impacts the image reconstruction of black holes [76]. Arguably the most serious systematic that must then be confronted in LCS studies is how to reliably extract a targeted distribution, while minimizing numerical artifacts and bias. Numerous sophisticated methods, such as the Backus-Gilbert [77], maximum entropy [78], and Bayesian reconstruction methods have been explored as tools to aid in PDF extractions [57,79] and other observables more generally [80,81].
A common heuristic to regularize the inverse problem at hand, both in the global fitting of inclusive scattering data [82][83][84][85][86][87] and analogous lattice calculations, is to supply additional information in the form of physically motivated PDF parameterizations. In particular, a parametric form can incorporate the known divergence/convergence of PDFs at small-/large-x and enforce any parton sum rules explicitly. We note parametric frameworks in the literature generally differ in what functional form is used to smoothly connect the two limiting x-space regimes. Inspired by the phenomenological forms of the global fitting community and to establish a benchmark for the alternate extraction methods that follow, we opt to first regularize the inverse problem by parameterizing the valence/plus quark distributions according to where P (x) is a smooth interpolating polynomial and N −1 v = B (α + 1, β + 1)+ k λ k B α+ 1 + k+1 2 , β + 1 ensures the valence quark sum rule 1 0 dx f qv/N x, µ 2 = 1 is satisfied; the normalization of f q + /N x, µ 2 is not fixed by a sum rule, and is left to float in our fits. Given the limited range of Ioffe-time in our results, we will find the simplest 2-parameter ansatz with P (x) = 1 cannot be avoided. Where possible, the bias introduced by this   [82] and JAM20 [88], and the NNLO analyses of MSTW [89] and NNPDF [87] at the same scale.
highly-constraining choice will be studied by supplementing P (x) with additional halfinteger powers of x: P (x) = 1 + k λ k x (k+1)/2 , thereby increasing the flexibility of our parameterizations beyond the nominal PDF behavior We start with two-and three-parameter PDF parameterizations, where in the latter we take P (x) = 1 + δx. The cosine/sine transforms of the PDF forms (4.7) and (4.8) are fit to the real/imaginary ITD data using first an uncorrelated least-squares regression with σ 2 Q the variance of the ITD, and {ν min , ν max } representing potential cuts on the data. These uncorrelated fits include all z/a ∈ {1, · · · , 12} and ap z ∈ {1, · · · , 6} × 2π/L.   [82] and JAM20 [88], and the NNLO analyses of MSTW [89] and NNPDF [87] at the same scale.
by the CJ [82] and JAM [88] collaborations, and three flavor NNLO determinations of MSTW [89] and NNPDF [87]. 1 Apart from the z/a ≥ 9 data, such an uncorrelated fit would seem to well describe the Re Q ν, µ 2 data and lead to valence PDFs that feature many structural similarities with phenomenological determinations at the same scale. The statistically consistent figure of merit for the two-and three-parameter fits tabulated in table 2, however indicates the data cannot distinguish between these models. The two-parameter fit to Im Q ν, µ 2 is clearly more heavily constrained by the z/a 7 data, and all but avoids points of the ITD originating from larger separations. Above x ∼ 0.4 the extracted f q + /N x, µ 2 n=2 C result likewise shares structural similarities with the shown phenomenological results. The lack of any large-ν constraint provided by Im Q ν, µ 2 entails a generally unconstrained fitted PDF in the small-x regime, although this relation is not bijective.

Direct extraction of PDFs from reduced pseudo-ITDs
A separate, though in principle equivalent, route to extract PDFs from these data is to directly apply the factorized relationship (2.5) having substituted the definition of the ITD (4.6): By assuming a PDF parameterization and performing a maximum likelihood regression of the double convolution and M ν, z 2 , the introduction of additional systematic errors from the evolution/matching steps and a potentially incorrect functional description of M ν, z 2 when interpolating its ν-dependence (e.g. eq. (4.5)) can all be avoided. The direct matching relationship between the PDFs and the reduced pseudo-ITD is given by where the one-loop kernels that match the {x, µ 2 }-dependencies of the valence/plus quark PDFs to the reduced pseudo-ITD are given by with the Altarelli-Parisi and scheme matching kernels modified tõ and Si (y) /Ci (y) are the integral sine/cosine functions and 3 F 3 (111; 222; −iy) is a generalized hypergeometric function [47,90]. A notable challenge of this direct approach is that accurate computation of the generalized hypergeometric function requires multi-precision arithmetic, which is computationally inefficient and invariably slows the parameter optimization. PDFs extracted by directly performing parametric fits (4.13) to the reduced pseudo-ITD data are denoted type-K.

JHEP11(2021)148
(a) (b) Figure 10. Two-parameter valence (a) and plus (b) quark PDF resulting from type-C (red) and type-K (indigo) fits to the unpolarized nucleon ITD and reduced pseudo-ITD, respectively. The direct matching fits are consistent with the cosine/sine transform of the model PDF fit to the ITD.  Simple two-parameter PDFs obtained from uncorrelated type-K fits are shown in figure 10a and figure 10b, together with the same phenomenological determinations and the uncorrelated type-C two-parameter PDF fits. The type-C and type-K fits are statistically consistent. This is confirmed by comparing the type-K fit results in table 3 to the type-C results in table 2. However, the central values of the type-K fits suggest that at small-x the f qv/N (x) is more divergent and the f q + /N (x) is instead convergent for small-x at the scale µ = 2 GeV. The factor of two or three increase in the figure of merit when switching from type-C to type-K fits is the first indication of puzzling behavior in M ν, z 2 . We reiterate the naive two-parameter fits capture the known limiting regimes of the PDFs. The poor figures of merit in the type-K fits hint that M ν, z 2 at this stage apparently does not align well with expectations from the direct matching (4.13). This is a potentially disastrous conclusion. To gain some insight, we now consider the data correlations.

Data correlation
The data featured in this work, and indeed any lattice calculation, naturally are correlated. By beginning this subsection with uncorrelated fits, we highlight that without knowledge or through the simple neglect of data correlations, which appears to be common in the literature, one might incorrectly assume an adequate description of the data has been JHEP11(2021)148 Focus is given to the small-ν region. The correlated two-parameter fit is seen to deviate appreciably from the precise small-ν data.
achieved. These correlations must be taken into account in order to provide a rigorous accounting of mutual fluctuations in the data and thus an agnostic PDF determination.
Simply repeating the two-parameter fit to Re Q ν, µ 2 , only this time accounting for the data covariance Cov with q k = Q ν, µ 2 k − Q fit ν, µ 2 k , we arrive at a much different conclusion shown in figure 11a. The visual discrepancy between the ITD and two-parameter fit is stark, and leads to a correlated figure of merit of O (40). Although the fit misses nearly all of the moderate to large-z points, figure 11b illustrates the large increase in the figure of merit is primarily due to the slight deviation from the very precise z/a 4 data.
Visualizing the data covariance in the real component in figure 12a, it is clear the low-momentum data ap z ≤ 4π/L are strongly correlated amongst each other and correlate weakly with the ap z ≥ 8π/L data; some mild correlation is visible with the ap z = 6π/L data with z/a ≤ 6. Within the ap z ≤ 4π/L channels the strongest correlation can be found in the shortest Wilson line data. These observations provide an explanation for the poor correlated two-parameter fit in figure 11a -the strongest correlation is with the most precise data in our calculation causing any correlated fit to favor the small-ν data. Indeed strong correlation is also observed amongst the momentum channels ap z = {4, 5, 6} × 2π/L, but the signal-to-noise degradation for these high-momentum data minimizes their effect on any fit. It is interesting this delineation corresponds to the transition from an unphased to phased eigenvector basis. The data covariance in the imaginary component, shown in figure  adjacent Wilson line lengths (e.g. z/a = 4 and z/a = 5). It is then no surprise that a correlated two-parameter PDF parameterization of Im Q ν, µ 2 is also met with a poor figure of merit. The non-trivial structures of correlation evident in these data are indicative of our simple PDF parameterizations (4.7) and (4.8) being inappropriate for these data. The above puzzling, and indeed worrisome, conclusions are given a deeper quantitative understanding in the following sections.

Classical orthogonal polynomials
The phenomenological parameterizations we have considered thus far are but one way to regulate the ill-posed inverse relation between the ITD/reduced pseudo-ITD and the corresponding PDF. These parameterizations nevertheless introduce a model dependence into the extracted PDF. Any PDF faithfully reported from a lattice calculation should take into account the space of functions that smoothly connects the x → 0 and x → 1 limits. As an alternative means to describe the valence/plus quark sectors and minimize model bias, we propose to parameterize the PDFs by a complete basis of classical orthogonal polynomials [55]. The leveraging of orthogonal polynomials to obtain an unknown distribution is not unique to this work. The approach we adopt parallels efforts to extract PDFs from phenomenological fits of inclusive processes [84,91], as well as distribution amplitudes [92][93][94] and inelastic scattering cross sections [95] from matrix elements calculated in lattice QCD. Consider the Jacobi (hypergeometric) polynomials As the set of polynomials {Ω (α,β) n } span x ∈ [0, 1], a PDF can be expressed generically as with expansion coefficients C (α,β) a,n . The parameters {α, β} lose their familiar characterization of the x → 0/x → 1 PDF behaviors in place of delineating between different choices of bases. The expansion in Jacobi polynomials in eq. (4.19) is thus entirely generic and model-independent. However, the series of Jacobi polynomials must in practice be truncated at some finite order, n. The bias then introduced may be studied by fixing the order of truncation and determining the optimum {α, β}, or tuning {α, β} to capture generic properties of a PDF and subsequently optimize the order of truncation -we will adopt the former.
Our strategy to parameterize the reduced pseudo-ITD using a set of Ω (α,β) n (x) will be met by similar numerical difficulties as the type-K fits discussed above. The numerical effort is lessened by considering a Taylor series expansion in ν for fixed separations z 2 of the direct matching kernels K v/+ xν, z 2 µ 2 and eq. (4.19). The contribution of an n th -order Jacobi polynomial Ω (α,β) n (x) to Re M ν, z 2 and Im M ν, z 2 is given by Expanding the direct matching kernels K v/+ xν, z 2 µ 2 in even/odd powers of ν one finds

JHEP11(2021)148
and the constants γ n and d n are the leading order moments of the Altarelli-Parisi and scheme matching kernels derived in [96], respectively. The sum over k is to be performed to assure convergence for a given value of ν -we have identified k max = 75 as providing more than adequate numerical precision, in reasonable computation time. With the above definitions, the leading-twist valence and plus quark PDFs describe the reduced pseudo-ITD components according to where the C lt (α,β) v/+,n are the Jacobi polynomial expansion coefficients. Since the space-like matrix element (2.1) is on-shell and given our choice of an O (a) improved lattice QCD formulation, any contaminating effects of this calculation must scale as a 2 Λ 2 QCD , z 2 Λ 2 QCD , and as some ratio of the lattice spacing a and Wilson line length z. Apart from the Ioffe-time ν = p · z, these terms are the only possible dimensionless combinations of dimensionful parameters in this calculation, embodying discretization errors of the reduced pseudo-ITD that vanish in the continuum limit and higher-twist errors that survive the continuum limit. Since this calculation has been performed at only a single lattice spacing, the a 2 Λ 2 QCD effect cannot be quantified, leaving the nominal polynomial corrections to the short-distance factorization of the reduced pseudo-ITD (2.5) and the short-distance error arising from the bilocal operatorO  Rather than make a selection for the form of a potential discretization effect, we allow our data to dictate the precise form. This shall be discussed in the ensuing paragraphs. As the Fourier transform in ν of the reduced pseudo-ITD only has support on the momentum fraction interval x ∈ [−1, 1] [36], the contaminating effects we have discussed must also only have support in this interval. By construction, these corrections must be functions of ν 2 and ν in the real and imaginary components of M ν, z 2 , respectively. The same basis of Jacobi polynomials in (4.23) and (4.24) can then be used to parameterize the Ioffe-time dependence of these contaminating effects, with x-space distributions of the same form as (4.19) but with distinct expansion coefficients. We denote the coefficients of the corrections as C corr (α,β) τ,n , with τ = {v, +} indicating whether the effect arises in the valence/plus quark PDFs. Supposing, for simplicity, these effects enter at tree-level, their contributions to the reduced pseudo-ITD signals with σ  where κ corr is the dimensionless parameter which describes the scaling of each correction (e.g. κ corr = z 2 Λ 2 QCD ). Visualizing the Taylor-expanded matching kernels (4.20) and (4.21) at tree-level across a range of Ioffe-times in figure 13a and figure 13b, it is seen the polynomials σ (α,β) 0,n , η (α,β) 0,n reach an extremum in Ioffe-time commensurate with the polynomial order and asymptotically approach zero. This conveniently reflects the correct large-ν behavior of the ITD in the same limit (see eq. (4.4)). Since M 0, z 2 = 1 by construction, all corrections must vanish at zero Ioffe-time. Of the Jacobi polynomials used to describe the corrections at tree-level, only σ (α,β) 0,0 (0) = 0 (blue curve of figure 13a). Hence the corrections to the real component of M ν, z 2 shown in eq. (4.26) are restricted to order n ≥ 1.
The complete functional forms we apply to each component of M ν, z 2 are:

JHEP11(2021)148
The leading-twist (C lt τ,n ) and discretization (C az τ,n ) corrections are accompanied by twist-4 (C t4 τ,n ) and twist-6 (C t6 τ,n ) corrections, where we note each higher-twist correction is a hypercubic invariant. The twist-6 corrections have been included for exploration purposes despite the likelihood of their effect being minimal in the range of Ioffe-time for which we have statistically clean data (ν ∼ 10). We point out the discretization effect in (4.28) we adopt herein scales as a/ |z|; the correlated figure-of-merit of this fit to Re M ν, z 2 was found to be several times smaller when a/ |z| was used in place of a 2 /z 2 . In this sense, the reduced pseudo-ITD data was found to favor a short-distance discretization effect of the form a/ |z| in both real and imaginary components. Indeed other correlations of the lattice spacing with z (e.g. a 2 /z 2 ) may be present in our data, but are evidently sub-leading to a/ |z| which we consider for the remainder of this work.
There are further sources of systematic error in M ν, z 2 that have been neglected by our parameterizations (4.28) and (4.29). Were the full NLO coordinate space matching kernel considered in the expansion of each correction in Jacobi polynomials (cf. eq. (4.26) & (4.27)), it becomes clear the discretization and higher-twist corrections we consider would enter at NLO with a more complicated dependence on z, namely an additional logarithmic dependence on z 2 . These additional terms, however, each receive a further suppression by a factor of O (α s /π) ∼ 0.1, and are on the order of a few percent in the range of ν for which we have computed M ν, z 2 . Only for considerably large-z will the O (α s ) components of 0,n , become relevant. These additional sources of z 2 contamination are avoided, and hence excluded from (4.28) and (4.29), by noting the short-distance factorization (2.5) would be invalidated by such large separations.

PDF results with Jacobi polynomials
In the fits we perform according to (4.28) and (4.29), we elect to fix the order of truncation for the leading-twist and each type of correction, and numerically search for the optimal {α, β, C corr τ,n }. Treating each fitted parameter as non-linear in a maximum likelihood fit leads to wildly unstable results. The way forward is to recognize α, β are the only fitted parameters that are truly non-linear; the correction coefficients C corr τ,n are all linear. A maximum likelihood fit of the posterior distribution of the linear terms is then Gaussian and cheap to obtain.
The reduced pseudo-ITD fits of eq. (4.28) and eq. (4.29) are implemented with the help of the Variable Projection (VarPro) algorithm [97] for separable non-linear optimization problems. This reduces the dimension of the non-linear optimization from d = N corr + 2 to d = 2, where N corr are the number of linear correction coefficients. In our case, minimization is performed in the d = 2 Jacobi polynomial basis {α, β}, and any correction terms C corr (α,β) τ,n are solved for exactly in terms of the non-linear basis functions {σ We note without VarPro, optimizations with N corr ≥ 4 are numerically unstable, regardless of the type of correction included in the model.
Care needs to be taken as correction terms are included in eqs. (4.28) and (4.29), as physical insight can quickly be replaced with over fitting. The first sensible restriction to impose is for all x-space corrections O (a/ |z|), O(z 2 Λ 2 QCD ) and O(z 4 Λ 4 QCD ) to be sub-leading relative to the leading-twist PDF. It would be alarming to obtain, say, a twist-4 contribution JHEP11(2021)148 that is larger than the leading-twist PDF, given that the short-distance factorization of the pseudo-distributions implies leading-twist dominance. Such disastrous scenarios are avoided with several Bayesian constraints of a Gaussian form. So as to allow the reduced pseudo-ITD to dictate the best fit results, all Bayesian priors are fixed to zero. The hierarchy we desire is realized with the following prior widths: • Leading-twist: The validity of the entire Jacobi polynomial parameterization is guaranteed using shifted log-normally distributed priors to ensure α, β > −1. In practice, the log-normal prior on beta is shifted to β = 0 to secure β > 0 and hence convergent PDFs at x = 1.
As leading-twist and correction terms are added, the question becomes at which order each series of Jacobi polynomials should be truncated. We address this by scanning over all possible combinations of truncation orders for n lt ∈ {3, 4, 5, 6}, n az ∈ Z 4 , n t4 ∈ Z 5 , n t6 ∈ Z 3 , where n * are the orders of truncation in the fits (4.28) and (4.29). Figure 14 illustrates, for a rather large number of Jacobi polynomials {n lt , n az , n t4 , n t6 } = {6, 3, 4, 2}, the covariances of α, β and each linear correction term C corr (α,β) τ,n in fits to the real (4.28)

JHEP11(2021)148
and imaginary (4.29) reduced pseudo-ITD components for Wilson line lengths z/a ≤ 12. The covariance of each pair of fitted parameters is estimated via jackknife resampling (4.30) where fit parameters associated with jackknife sample n are denoted by f n,k , with jackknife averagef k . Without observing the quality of agreement between each fit and the reduced pseudo-ITD, it is clear several parameters correlate weakly or not at all with other parameters in the fit. This implies these weakly correlated parameters are not well-constrained by the data, and their removal will not affect the information content of the fit. For instance, the real component fit parameter covariances, shown in figure 14a, suggest the leading-twist expansion coefficients C lt n lt are constrained by the data for n lt ≤ 3, while C lt 4 , C lt 5 weakly correlate with the remaining parameters. The discretization, twist-4 and twist-6 corrections exhibit mild correlation for C az 1 , C t4 2 , C t4 3 , C t6 1 , with the remaining correction parameters largely unconstrained. In the imaginary component, the fit parameter covariances shown in figure 14b suggest a more nuanced pattern of correlation. Several leading-twist Jacobi polynomials appear to be well-constrained by the data, while the relative correlation between the C az and C t4 , C t6 parameters is increased relative to the corresponding entries in figure 14a.
This exercise demonstrates an important point. Although the VarPro implementation of the Jacobi polynomial fits allows for arbitrarily many leading-twist and correction coefficients, the reduced pseudo-ITD data simply do not contain enough information to constrain so many parameters. One should then expect the likelihood function is maximized for the real component of the reduced pseudo-ITD with truncation orders n lt ∼ 3 and n az ∼ 1 − 2, and n lt ∼ 3 and n az ∼ 3 for the imaginary component.
By scanning over the order of truncation for the leading-twist and correction terms parameterized by Jacobi polynomials, we find the likelihood of the functional (4.28) to describe Re M ν, z 2 with z/a ≤ 12 to be maximized for {n lt , n az , n t4 , n t6 } v = {4, 1, 3, 2} v . Likewise, the likelihood of the functional (4.29) to describe Im M ν, z 2 with z/a ≤ 12 is maximized for {n lt , n az , n t4 , n t6 } + = {3, 3, 1, 0} + . The fit results for each and their respective figures-of-merit are given in table 4. The Jacobi polynomial fits of orders {n lt , n az , n t4 , n t6 } v = {4, 1, 3, 2} v and {n lt , n az , n t4 , n t6 } + = {3, 3, 1, 0} + applied to the real and imaginary components of M ν, z 2 are presented in figure 15 and figure 16, respectively.
Considering first the real component fit, each set of Re M ν, z 2 for z/a ≤ 8 are well represented by the expansion in Jacobi polynomials. The main exception is the highest momentum point p z = 6 × (2π/aL) ∼ 2.47 GeV. The Re M ν, z 2 data for z/a > 8 are also reasonably well described, however the highest two momenta are seen to deviate. This behavior is not surprising despite the twist-4 and twist-6 corrections, which capture large-z 2 deviations, as the highest momentum data are subject to loss of signal in both the twoand three-point functions. The associated fit parameter covariances shown in figure 17a demonstrate the leading-twist, discretization and twist-4 corrections are well constrained by the Re M ν, z 2 data; as expected, the twist-6 corrections are only weakly constrained.    Table 4. Various Jacobi polynomial fits to the real and imaginary components of the unpolarized reduced pseudo-ITD for z/a ≤ 12. Each column represents distinct orders of truncation in the Jacobi polynomial expansions to the leading-twist, discretization, twist-4 and twist-6 corrections. The real and imaginary component fits were found to have the highest likelihoods of describing the data with truncation orders {4, 1, 3, 2} v and {3, 3, 1, 0} + , respectively. The dramatic effect even a single discretization term has on each fit is shown in the columns {4, 0, 3, 2} v and {3, 0, 1, 0} + .
The resultant leading-twist PDF f qv/N (x) and x-space distributions corresponding to the {n lt , n az , n t4 , n t6 } v = {4, 1, 3, 2} v Jacobi polynomial fit are gathered in figure 18b. As expected, the corrections in x-space are sub-leading to the leading-twist PDF. The Jacobi-parameterized leading-twist PDF, however, features many structural differences with the included phenomenological PDFs and the uncorrelated two-parameter PDF fit. Most evident is the softer approach to x = 1. Due to the valence quark sum rule, this enhances the low-to moderate-x region and leads to further tension with the phenomenological results. By evaluating the cosine transform of the pure leading-twist component, we see in figure 18a that the z/a 7 ITD data deviate successively further from the derived leading-twist ITD shown in purple. Whereas the uncorrelated two-parameter ITD fit shown in red attempts to capture all the z/a ≤ 12 data and indeed the unwanted impact from higher-twist effects, the Jacobi polynomial parameterization has effectively isolated and removed these polynomial-z 2 effects, leaving the pure leading-twist contribution. The quality of the Jacobi polynomial fit to the imaginary component of the reduced pseudo-ITD shown in figure 16 is more puzzling. The z/a ≤ 4 appear reasonably well represented by the expansion in Jacobi polynomials, but by z/a = 5 it is evident the data for a given z 2 segregate into two distinct groups -one for lattice momenta p latt ∈ {1, 2, 3} and another for p latt ∈ {4, 5, 6}. This distinction coincides with the switch from an unphased eigenvector basis to the phased bases ζ ± defined in eq. (3.7). The fit parameter covariances shown in figure 17b demonstrate a milder constraint of the first and second order leadingtwist Jacobi polynomials compared to the best fit of the real component. The discretization and twist-4 corrections are also seen to be well constrained by the data. The resultant leading-twist plus quark PDF f q + /N (x) and x-space distributions corresponding to the {n lt , n az , n t4 , n t6 } + = {3, 3, 1, 0} + Jacobi polynomial fit are illustrated in figure 19b. As in the real component fit, the corrections are sub-leading to the leading-twist PDF, which in this case is in agreement with the NNPDF result [87] for x ≥ 0.5. At small values of x, the leading-twist PDF parameterized by Jacobi polynomials is generally consistent with the two-parameter uncorrelated PDF fit. The sine transform of the pure leadingtwist component is shown in figure 19a together with the Im Q ν, µ 2 data at 2 GeV. Unlike the real component of the derived leading-twist ITD in figure 18a, the derived imaginary component of the leading-twist ITD does not agree with the Im Q ν, µ 2 data for any of the ap z 4π/L data with z/a 7. As the imaginary component of M ν, z 2 is optimally fit with three discretization corrections and only one higher-twist term, the {n lt , n az , n t4 , n t6 } + = {3, 3, 1, 0} + fit would suggest the imaginary component of the ITD is susceptible to less higher-twist effects in exchange for greater discretization effects. This is a tenuous conclusion, however, in light of the segregation of the Im M ν, z 2 data into two distinct clusters, a low-and high-momentum set, for large Wilson line lengths. A future study exploring the side effects of phased distillation is warranted.
By far the biggest indicator of a reasonable description of the reduced pseudo-ITD data via Jacobi polynomials is a discretization term. Repeating the above Jacobi polynomial fits but leaving out any discretization corrections, namely {n lt , n az , n t4 , n t6 } v = {4, 0, 3, 2} v and {n lt , n az , n t4 , n t6 } + = {3, 0, 1, 0} + , the correlated figures of merit increase considerably to unacceptable values (see table 4). This same conclusion is reached when cuts on momentum and Wilson line lengths are made. Since the discretization term we have included is of O (a/ |z|), its effect is most pronounced at short distances. This is precisely the regime JHEP11(2021)148 The distributions are compared with the uncorrelated 2-parameter phenomenological fit of eq. (4.8) (red), as well as the NLO global analyses of CJ15 [82] and JAM20 [88], and the NNLO analyses of MSTW [89] and NNPDF [87] at the same scale.
wherein the short distance factorization (4.1), or equivalently (4.12), is applicable. This motivates a more detailed look at the short-distance behavior of the computed reduced pseudo-ITD.

On the numerical consistency with DGLAP
The one-loop matching relationship between the ITD and the reduced pseudo-ITD (4.1) implies that Q ν, µ 2 = M ν, z 2 at tree-level. The scatter that exists for a given z 2 should ideally be compensated at O (α s ) by the ln z 2 -dependence produced by the DGLAP evolution, up to large-z 2 higher-twist corrections. In this section we study the z 2 -dependence of M ν, z 2 more closely, and investigate whether the observed dependence is numerically consistent with DGLAP, thus yielding a truly z 2 -independent ITD. We begin by focusing on the real component of the reduced pseudo-ITD. The dependence of Re M ν, z 2 on the invariant space-like interval z 2 can be most easily visualized by parameterizing the valence pseudo-PDF P v x, z 2 ; α, β by a simple two-parameter phenomenological form and fitting its cosine-transform to Re M ν, z 2 separately for each z 2 . In order to more readily expose the role of the Altarelli-Parisi kernel, we impose the added restriction β = 3. This choice not only captures the naive x → 1 behavior of the nucleon's valence quark PDF [98], but also forces α to reflect any z-dependence in the reduced distribution; further, this value of β is in statistical agreement with those obtained from the uncorrelated ITD fits (see table 2). Figure 20 illustrates the cosine-transform of the model valence pseudo-PDF (5.1) fit separately to each z 2 of the real component of the reduced pseudo-ITD. The cosinetransforms of P v x, z 2 ; α, β = 3 are seen to describe Re M ν, z 2 quite well for z/a 10, with the greatest tension seen for the highest momentum point for each separation. The fits for z/a ≥ 13 are also shown for completeness, but are clearly noise dominated. Also noteworthy, the highest figures-of-merit are observed for the smallest separations, with a somewhat monotonic reduction until z/a 11. The dependence of the fitted values of α on the separation z/a is visualized for Re M ν, z 2 in figure 21. As a function of z/a, α decreases with the Wilson line length, matching expectations from the Altarelli-Parisi evolution of the pseudo-PDF. However, it is clear Re M ν, z 2 depends linearly on z/a for z/a 12, most notably for small-z. This manifest lack of ln z 2 behavior of Re M ν, z 2 at short distances immediately suggests tension with the presumed DGLAP evolution of the pseudo-PDF. To determine if this z 2 -dependence in Re M ν, z 2 is nevertheless numerically consistent with DGLAP, the one-loop matching relationship between the reduced pseudo-ITD and ITD is applied. In the ideal scenario where the z 2 -dependence of Re M ν, z 2 is exactly described by DGLAP, the matched ITD will be independent of the interval z 2 up to polynomial corrections for large-z 2 . Rather than perform the matching step to a common scale in MS using a smooth polynomial in Ioffe-time (e.g. eq. (4.5)) as was done in section 4, we leverage the cosine-transform of the model pseudo-PDF (5.1) as the smooth and continuous description of the reduced pseudo-ITD data. That is, we perform the matching of Re M ν, z 2 to a common scale in MS according to 2) where P uν, z 2 ; α, β = 3 is the cosine-transform of the model pseudo-PDF P v x, z 2 ; α, β = 3 expressed in a closed form by a generalized hypergeometric function For an explicit, albeit crude, conversion to MS, we set α = 0.2 in eq. (5.3). Our strategy to expose any z-dependence in the ITD Q ν, µ 2 remains identical to the reduced distribution above. The resultant matched ITD at 2 GeV in MS is once more fit using the two-parameter form in eq. but rather the valence PDF itself. As illustrated in figure 22, each z 2 of the matched ITD is well described by the simple two-parameter form. The poorest figures-of-merit are again observed for the smallest (z/a 3) and largest (z/a 13) separations. The dependence of the fitted values of α on the separation z/a for Re Q ν, µ 2 is illustrated in figure 23. For 4 z/a 11 the fitted value of α is observed to be independent of z/a and hence numerically consistent with DGLAP in said interval. Remarkably, however, the values of α for the shortest separations, namely z/a 4, deviate increasingly from this constancy as z/a → 1. A subsequent analysis of the imaginary component of both the reduced pseudo-ITD and matched ITD arrived at a similar conclusion, but has been omitted for brevity.

Jacobi polynomial corrections -discretization effects
The findings above rigorously demonstrate the reduced pseudo-ITD is numerically inconsistent with DGLAP in the small-z regime. Whether matching the reduced pseudo-ITD to the light-cone ITD or directly to the light-cone PDF, the presence of the Altarelli-Parisi denote Re M ν, z 2 ≡ Re M ν, z 2 − Re M az ν, z 2 . Based on figure 24a, the removal should shift the small-z points of Re M ν, z 2 to larger values, with the largest impact for ν ∼ 4.5. Figure 25 juxtaposes the original Re M ν, z 2 and discretization corrected Re M ν, z 2 in the interval ν ∈ [0, 2.5]. Although the differences are numerically small, at small Ioffe-times Re M ν, z 2 is noticeably larger than the uncorrected reduced pseudo-ITD. The importance of removing this discretization effect is quantitatively discerned by repeating the DGLAP investigation for Re M ν, z 2 .
Parameterizing the discretization corrected valence pseudo-PDF P v x, z 2 with the two-parameter form in eq. (5.1) and fitting its cosine-transform to Re M ν, z 2 with β = 3, the z-dependence of Re M ν, z 2 is once more reflected in the variation of α with z/a. As illustrated in figure 26a, α now varies non-linearly with z/a for z/a 4 and linearly for 4 z/a 11. Whether this markedly distinct z-dependence (cf. figure 21) is numerically consistent with DGLAP is once again checked by performing the matching to a common scale in MS using eq. (5.3), and repeating the two-parameter fits to the discretization corrected ITD Re Q ν, µ 2 for each z 2 and with β = 3. The resulting fitted values of α are presented in figure 26b. Relative to the z-dependence of the uncorrected ITD shown in figure 23, the variation of α with z/a is considerably more constant for z/a 11. In other words, the ITD is seen to fall into better agreement with DGLAP in the short-distance regime following removal of the O (a/ |z|) effect. That the optimal Jacobi polynomial fits {n lt , n az , n t4 , n t6 } v = {4, 1, 3, 2} v and {n lt , n az , n t4 , n t6 } + = {3, 3, 1, 0} + provide the best description of M ν, z 2 can now be quantitatively explained by the compensating effect the O (a/ |z|) term provides. The poor quality of the correlated phenomenological fits to the matched ITD, as well as the correlated Jacobi polynomial fits to M ν, z 2 without any corrections, are a direct result of attempting to fit a singular function in z to data that do not exhibit singular behavior. By excluding z/a 4 and z/a 11, the short-distance tension and any large-z polynomial effects can be removed yielding reduced pseudo-ITD or matched ITD data that are well in line with theoretical expectations. Such cuts are common in the literature, however their nominal effect is to neglect deviating behavior.
Although the DGLAP investigation has been shown for the real component of the reduced pseudo-ITD, the considerable reduction in the correlated figure-of-merit when discretization effects are included in fits to Im M ν, z 2 (e.g. {n lt , n az , n t4 , n t6 } + = {3, 3, 1, 0} + versus {n lt , n az , n t4 , n t6 } + = {3, 0, 1, 0} + in table 4) indicates the imaginary component of the raw reduced pseudo-ITD likewise deviates from expectations of DGLAP at short-distances. The central question left for future research is the origin of this discretization effect.

Conclusions
In this work we presented in detail the first lattice calculation of the unpolarized nucleon PDF in the distillation framework, showing its efficacy and emphasizing in detail the advantages of this smearing technique. We have observed that distillation can yield a pseudo-distribution of superior quality, addressing one of the principal components limiting the competitiveness of contemporary lattice QCD calculations of PDFs with respect to precision phenomenological extractions. For the unpolarized PDF, the impact of this lattice calculation on global fits was found to be rather marginal. The quality of data, however, hints that a future calculation of a distribution less constrained by experiment may benefit from the use of distillation. By validating the distillation method in the unpolarized nucleon PDF, our study opens new avenues of synergy between lattice and phenomenology in the spirit of [11,99,100]. We performed a careful study of the correct extraction of the matrix elements that yield the pseudo-ITD and we analyzed in detail its real and imaginary components. This was followed by a discussion of correlations of lattice data and how the commonly used uncorrelated fits can yield erroneous results and can severely flaw the extraction of the PDFs. We regulate the ill-posed inverse problem that relates the ITD to the PDF with the use of Jacobi polynomials, and stress the necessity of removing discretization

JHEP11(2021)148
Frontera [101] that have contributed to the results in this paper. We acknowledge the facilities of the USQCD Collaboration used for this research in part, which are funded by the Office of Science of the U.S. Department of Energy. This work was performed in part using computing facilities at the College of William and Mary which were provided by contributions from the National Science Foundation (MRI grant PHY-1626177), and the

JHEP11(2021)148
The inversion cost present in lattice calculations that implement a standard spatial smearing kernel must solve for the Green's function represented by where D is the chosen discretization of the Dirac operator, and the N c color and N s spinor components are denoted by Latin and Greek indices. Such point-to-all propagators stand in contrast to the so-called solution vectors required of a distillation based calculation: where ξ (k) is an eigenvector of the gauge covariant Laplacian with k th largest eigenvalue. The separation of the elemental (interpolator) construction from the solution vector computation with distillation implies the total inversion cost per configuration needed to realize matrix elements sensitive to the Ioffe-time pseudo-distribution is given by where N eigs is the rank of the distillation space, N tsep the number of three-point source-sink separations, N ζ the number of phased eigenvector bases (including unphased), and N src the number of temporal source origins. In this calculation we have selected N tsep = 6, N eigs = 64, N ζ = 3, and N src = 4, for which N dist inv /cfg = 16128. We invite the reader to compare this inversion cost with that needed to isolate matrix elements when using sequential source methods (discussed, for example, in ref. [109]), which are common cost saving measures when exploiting standard spatial smearing kernels.
One aim of this body of work was to demonstrate the benefits distillation can provide in calculations of hadronic structure. An earlier calculation of the isovector unpolarized valence quark content of the nucleon using the pseudo-distribution formalism on the a094m358 ensemble used in this work was presented in ref. [58]; this previous calculation instead used Jacobi smearing. Figures 27 and 28 overlay the real and imaginary components, respectively, of M ν, z 2 reported in ref. [58] with our results first presented in figures 4a and 4b. Evidently the pseudo-ITD obtained with distillation is considerably more precise than the Jacobi smeared data of ref. [58], with improvement factors at times exceeding an order of magnitude.
The data from ref. [58] was obtained using the sequential operator technique, for which the number of inversions needed per configuration is given by 4) with N zsep the number of Wilson line lengths. With N zsep = 17 and N ζ = 5 in ref. [58], it follows N Seq.Op inv /cfg = 8640. In other words, by roughly doubling the number of inversions per configuration, a step quite feasible with contemporary computing resources, we have achieved with distillation a reduction in error far greater than the naive 1/ √ 2 reduction factor. Furthermore, we have relied on N cfg × N srcs = 349 × 4 = 1396 measurements, while N cfg × N srcs = 417 × 8 = 3336 measurements were used in ref. [58]. The reader is directed to JHEP11(2021)148  [58] using Jacobi smearing (crosses). (b) Focusing on the small-ν region, the improvements afforded by distillation are seen to be considerable. Data from ref. [58] has been shifted horizontally for legibility.  [58] using Jacobi smearing (crosses). (b) Focusing on the small-ν region, the improvements afforded by distillation are seen to be considerable. Data from ref. [58] has been shifted horizontally for legibility. ref. [109] for further comparisons between distillation and Jacobi smearing in the calculation of nucleon charges.
An alternate inversion strategy, but no less ubiquitous, is that of the sequential propagator approach. In this case the source and sink interpolating operators are held fixed, with the series of inversions beginning at the source and running through the sink. In this manner, matrix elements of a variety of operator insertions and momentum transfers can be computed without additional inversions. The number of inversions needed per configuration is given by where N tsep are the number of source-sink separations, N mom the number of momentum projections at the sink, N proj the number of spin projectors applied to the nucleon interpolating field, N s the number of open spinor indices after projection of the nucleon interpolator, and N op the number of distinct interpolators at the sink. As standard three-point functions enforce three-momentum conservation on an ensemble average by performing explicit momentum projections only at the operator insertion and sink interpolator, the sequential propagator cost (A.5) only scales linearly with N mom and N op rather than also with the number of interpolators and momentum projections at the source. The factor of two in eq. (A.5) designates one sequential propagator for each of the two quark flavors in the nucleon. Translating these factors into our calculation with distillation using a single unpolarized non-relativistic nucleon interpolator, in the sequential propagator framework one would have N op = 1, N proj = 1 and N s = 2 for a total of N Seq.Inv.
inv /cfg = 4056. Indeed the sequential propagator strategy is often favored in place of the sequential operator framework discussed above and featured in ref. [58], as the inversion cost is roughly half and should produce identical correlation functions on an ensemble average. Implicit in the sequential operator cost (A.4) is the added cost per operator insertion and momentum transfer: 1 + N zsep N ins N q which reduces to (1 + N zsep ) in the case of the unpolarized PDF. The fixed inversion overhead for distillation represented by eq. (A.3) means the unpolarized, helicity and transversity PDFs are simultaneously accessible with the same set of inversions, while the inversion cost with the sequential operator method would increase by a factor of 1+3Nzsep 1+Nzsep . Conventional smearing and the associated sequential inversion methods are then far outpaced by distillation when one is interested in off-forward observables -the sequential operator approach would require an additional inversions from the unpolarized PDF case, while the sequential propagator method would require a prohibitive multiplicative increase by a factor N mom for the number of distinct sink momenta. In an analogous fashion, were one interested in a variational analysis of a matrix of three point correlators, the sequential propagator method would require a further multiplicative increase by a factor N ops , while the cost increase with the sequential operator method would manifest in additional contractions rather than inversions. This paradigm bears considerable weight on the few calculations of Generalized Parton Distributions (GPDs) from lattice QCD (e.g. [110][111][112]) present in the literature, wherein a large set of momenta and conceivably operator bases are required to map out the multi-dimensional structure of GPDs while effectively controlling excited-states. The cost of sequential inversion methods used in tandem with conventional smearing clearly grows rapidly in the off-forward regime, while the inversion cost with distillation remains fixed at (A.3).
The numerical benefits presented herein and the cost saving prospects possible with distillation, especially in GPD calculations, strongly motivates its continued use in future calculations of hadronic structure.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.