Lattice continuum-limit study of nucleon quasi-PDFs

The quasi-PDF approach provides a path to computing parton distribution functions (PDFs) using lattice QCD. This approach requires matrix elements of a power-divergent operator in a nucleon at high momentum and one generically expects discretization effects starting at first order in the lattice spacing $a$. Therefore, it is important to demonstrate that the continuum limit can be reliably taken and to understand the size and shape of lattice artifacts. In this work, we report a calculation of isovector unpolarized and helicity PDFs using lattice ensembles with $N_f=2+1+1$ Wilson twisted mass fermions, a pion mass of approximately 370 MeV, and three different lattice spacings. Our results show a significant dependence on $a$, and the continuum extrapolation produces a better agreement with phenomenology. The latter is particularly true for the antiquark distribution at small momentum fraction $x$, where the extrapolation changes its sign.


I. INTRODUCTION
The calculation of parton distribution functions (PDFs) using lattice QCD has seen renewed interest in recent years [1][2][3][4], driven in part by the introduction of the quasi-PDF method [5,6]. This method requires nucleon matrix elements of a nonlocal operator containing a Wilson line, which must be computed on the lattice. Previous calculations of quasi-PDFs and related observables using the same operator by ETMC are given in Refs. [7][8][9][10][11][12][13][14][15] and by other collaborations in Refs. .
The presence of a Wilson line in the nonlocal operator introduces a power divergence. This divergence must be exactly removed by the renormalization procedure so that a finite continuum limit can be obtained. Furthermore, in contrast to the case of local operators, the use of a lattice action with exact chiral symmetry or at maximal twist does not eliminate all discretization effects linear in the lattice spacing a [45][46][47]. This means that in a lattice setup where most observables have only O(a 2 ) lattice artifacts, quasi-PDFs can nevertheless have O(a) contributions. For both of these reasons, it is important to numerically study the approach to the continuum limit so that future calculations will be better equipped to control all sources of systematic uncertainty.
There exist some previous studies using more than one lattice spacing. Ref. [45] includes an early analysis using two of the three lattice spacings used in this work. Nonperturbative renormalization was studied using two lattice spacings in Ref. [30], and the same two lattice spacings were used for studying pion PDFs in Ref. [42]. After the first version of this paper was submitted, two more works appeared. Ref. [48] presents a study of nucleon PDFs using three lattice spacings and three different pion masses, in which the lowest two pion masses were each studied using a single lattice spacing and the highest pion mass was studied using two lattice spacings. Finally, zero-momentum pion matrix elements were computed in Ref. [49] using multiple actions and up to four lattice spacings per action.
In this paper, we present a study of the approach to the continuum limit of isovector nucleon unpolarized and helicity parton distributions using three lattice ensembles, each having a different lattice spacing but with otherwise similar parameters. Section II describes the ensembles and the observables we compute. A dedicated study on one ensemble of systematic effects from excited-state contamination is reported in Section III. Renormalization factors are arXiv:2011.00964v3 [hep-lat] 13 Apr 2021 Name β aµ l size a (fm) mπ (MeV) pzL/(2π) pz (GeV) ts/a ts (fm) N conf Nsamp A60 1.90 0.0060 24 3 × 48 0.0934 (13) I. Parameters of the three N f = 2 + 1 + 1 lattice ensembles: gauge coupling β, bare light quark mass aµ l , and size. The pion mass mπ and lattice spacing a (determined via the nucleon mass) are taken from Ref. [50]. Nucleon three-point functions are computed with momentum p = (0, 0, pz) and source-sink time separation ts. The total number of gauge configurations is given by N conf ; on each one, we use an evenly-spaced grid of 32 source positions, with a random overall displacement, yielding Nsamp = 32N conf samples.
obtained using two different methods in Section IV; in addition, we study a ratio of matrix elements that cancels the renormalization. In Section V, we take the continuum limit, both for position-space matrix elements and for PDFs. Finally, conclusions are given in Section VI.

II. LATTICE SETUP
We use three lattice ensembles that differ primarily in their lattice spacings a = 0.0644, 0.0820, and 0.0934 fm. These have dynamical degenerate up and down quarks with pion mass approximately 370 MeV and dynamical strange and charm quarks with near-physical masses, i.e. N f = 2 + 1 + 1. The gauge action is Iwasaki [51,52] and the fermions use Wilson twisted mass tuned to maximal twist. These ensembles were generated by ETMC [53]; parameters for the three used in this work are given in Table I. The ensemble with intermediate lattice spacing, B55, was previously used by some of us for studying quasi-PDFs in Refs. [7][8][9].
Isovector quasi-PDFs are obtained from nucleon matrix elements of the nonlocal operator where bold symbols denote Euclidean four-vectors, ψ is the doublet of light quarks, W is a Wilson line, τ 3 selects the isovector u − d flavor combination, and we have chosen to extend the operator in the third spatial direction. We employ five steps of stout smearing [54] in the definition of W . The operator's nucleon matrix elements can be written as where µ represents the scale at which O is renormalized. Taking the Fourier transform, we obtain the unpolarized and helicity quasi-PDFs,q These are related to physical PDFs through factorization, and a similar expression applies to the helicity case.
The details of our calculation are similar to Ref. [12], although we use nucleon momenta only in the +ẑ direction and do not improve statistics by averaging over equivalent directions. The proton interpolating operator is defined using Wuppertal-momentum-smeared quark fields [55,56], with the smearing performed using APE-smeared gauge links [57].

III. EXCITED-STATE EFFECTS
On ensemble A60, we performed a dedicated study of excited-state effects by varying the source-sink separation t s /a from 4 to 10. The nucleon effective energy on this ensemble is shown in Fig. 1; although momentum smearing yields ts/a N conf Nsrc Nsamp  {4, 5, 6, 7} 315 4  1260  8  315 8  2520  9  315 16  5040  10  1260 32 40320   TABLE II. Statistics used for excited-state study on ensemble  a good signal at moderate source-sink separations, the statistical uncertainty still grows rapidly at large separations. Therefore, we use much larger statistics for the larger separations, as given in Table II. Matrix elements are obtained from two-point and three-point correlation functions C 2pt (t s ) and C Γ,z 3pt (τ, t s ), where t s is the Euclidean time separation between the source and the sink and τ is the Euclidean time separation between the source and O Γ (z). We consider two estimators for the matrix element h Γ (z): where S Γ,z (t s ) ≡ a ts/a−1 τ /a=1 and ∆E is the energy gap to the lowest excited state.
Results are shown for the unpolarized and helicity matrix elements in Figs. 2 and 3. For both observables the excitedstate effects are similar. In the real part at small z, the dependence on t s is weak, especially for the unpolarized case where h γ0 (0) is a conserved charge. For z > 6a, h eff dips below zero at small t s and this negative part is substantially reduced when t s is increased. In the imaginary part, the negative peak around z = 3a is reduced in magnitude when t s is increased.
For most values of z, h ratio eff (z) with t s = 10a is consistent with the value for t s = 8a and 9a and also with h summ eff (z) for t s = 5a and 6a. Therefore we conclude that excited-state effects are reasonably under control using the ratio method with the largest time separation, and we choose to use similar separations for the two other ensembles. However, the analysis in the rest of this paper differs slightly from the excited-state study: instead of simply taking the midpoint τ = t s /2 in C Γ,z 3pt (τ, t s ), we average over several central values of τ to reduce the statistical uncertainty. The resulting bare matrix elements for all three ensembles are shown in Fig. 4.
One risk of studying excited-state effects using just one ensemble is that insufficiently controlled excited-state contributions on the other ensembles could be mistakenly interpreted as discretization effects. To reduce this possibility, Effective hγ 0 versus z: real part (top) and imaginary part (bottom). For each z, the ratio-midpoint results are shown using seven source-sink separations, increasing from left to right (red), and the summation-method results are shown using the shortest three source-sink separations (blue).  t s was chosen to be slightly larger on the two ensembles that lack an excited-states study. Furthermore, our findings in Section IV C, that accounting for the leading effect of small differences in p z improves the approach to the continuum, and in Section V, that the dependence on a is typically monotonic, are both consistent with discretization effects and not excited-state effects playing the dominant role in this study.

IV. RENORMALIZATION
Renormalization of the nonlocal operator O(z) was a stumbling block in rigorously calculating quasi-PDFs and was absent in the earliest lattice QCD calculations [7,8,16,17]. In contrast with local quark bilinears that diverge logarithmically, O(z) contains a Wilson line that introduces a power divergence. In order to obtain a continuum limit, it is essential that this divergence be removed exactly, meaning that lattice perturbation theory is inadequate. Nonperturbative renormalization prescriptions [9,19,45], introduced more than three years after the first lattice quasi-PDF calculations, are necessary.
We employ two different methods for nonperturbative renormalization, both of which involve imposing renormalization conditions on Green's functions evaluated on Landau-gauge-fixed lattices 1 . For this, we use the N f = 4 twisted mass ensembles from Ref. [59] listed in Table III. These have the same action and bare coupling as the ensembles used for computing nucleon matrix elements. However, because of the difficulty in reaching maximal twist with four degenerate light fermions, we instead average over pairs of ensembles with opposite PCAC masses. After renormalizing O Γ (z) in a nonperturbative intermediate scheme, perturbation theory is used to convert first to the MS scheme and then to a modified MS (MMS) scheme [12]. The latter cancels a log(z 2 ) divergence in the MS-renormalized matrix element at short distance and enables a matching between quasi-PDF and PDF that conserves charge.
The first method is the whole operator approach, where renormalization conditions are imposed independently on O Γ (z) for each z, producing a separate renormalization factor for each z. The procedure is very similar to methods commonly used for local quark bilinears, and the nonperturbative intermediate scheme is RI -MOM.
The second method is the auxiliary field approach, where the nonlocal operator is rewritten as a pair of local operators in an extended theory. Renormalization conditions are imposed on those local operators and on the action of the extended theory, producing a minimal set of renormalization parameters. The nonperturbative intermediate scheme, RI-xMOM, uses a mixture of momentum space and position space.
When z = 0, it is a special case where O Γ (z) is a local operator, namely a vector or axial current. For this point, we use the renormalization factor for the corresponding local operator determined in Ref. [59]. In the next two subsections we discuss each method and their sources of systematic uncertainty. In a third subsection, we form a ratio of nucleon matrix elements to cancel the renormalization of O Γ (z) and study the continuum limit of the ratio.

A. Whole operator approach and RI -MOM scheme
The Rome-Southampton approach [60] and its RI ( ) -MOM schemes are commonly used to determine renormalization factors of local operators. Our prescription for the nonlocal operator O Γ (z) closely follows Refs. [9,61] and the improvements from Ref. [12] for controlling systematic uncertainties; we refer the reader to those references for a more detailed discussion.
In Landau gauge and in momentum space, we compute the fermion propagator S q [Eq. (13)] and the amputated vertex function V O , with the operator O inserted at zero momentum transfer. We impose the conditions at each value of z, where X Born is the tree-level value of X. As a shorthand, we write Z V for the renormalization of the unpolarized operator O γ0 and Z A for the helicity operator O γ3γ5 . We choose the RI renormalization scale, µ 0 , so that the vertex momentum p has the same components in all spatial directions, that is, ap = 2π Ls (n t + 1 4 , n, n, n) with integer n and n t . More precisely, we choose momenta with P 4 ≡ ( i p 4 i )/( i p 2 i ) 2 ≤ 0.32, in order to suppress finite-a effects that break rotational symmetry [59,62].
The renormalization factors are calculated on the N f = 4 ensembles given in Table III; these have the same bare coupling β as the N f = 2 + 1 + 1 ensembles used for the bare matrix elements. The renormalization procedure can be summarized in the following steps:  2. Averaging of the two ensembles at opposite aµ sea PCAC values followed by chiral extrapolation of the form Z 0 + (aµ)Z 1 (or quadratic in am π ) for each lattice spacing. For β = 1.90 we take the average of the four ensembles, as there is only one aµ value available. For all three β values, we find a very mild dependence on the pion mass, similarly to what was found for other ensembles [12].
3. Conversion to the MS scheme and simultaneous evolution to the scale 2 GeV, using the expressions from Ref. [61]. 4. Elimination of residual dependence on the RI scale by fitting to extrapolate (aµ 0 ) 2 → 0. An extensive study on the choice of the renormalization scale and the corresponding systematic uncertainties can be found in Ref. [12]. The optimal fit range for all β values is (aµ 0 ) 2 ∈ [1, 3].
5. Conversion to the MMS scheme, which is necessary in order to apply a matching formula that satisfies particle number conservation.
The final estimates for renormalization factors are shown in Fig. 5. For the real part, the results with β = 1.90 and 1.95 are very similar, but the latter has a smaller imaginary part. The finest lattice spacing, β = 2.10, has a larger real part. The renormalized matrix elements from the three lattice spacings are shown in Fig. 11 and their approach to the continuum limit is discussed in Section V A.

B. Auxiliary field approach and RI-xMOM scheme
The auxiliary field approach [45,[63][64][65] introduces a new field ζ(z) whose propagator is a Wilson line along thê z direction. This allows the nonlocal operator in QCD to be represented using the local operator φ ≡ζψ in the extended theory: The problem becomes that of renormalizing the action for ζ and the composite operator φ; one finds that three parameters are sufficient to renormalize all operators O Γ (z) [45]: Filled blue squares, orange circles, and green diamonds show the data with stout-smeared links on the coarse, medium, and fine lattice spacings, respectively. Diamonds with black outlines show the data on the fine lattice spacing without smearing. Note that a hypercubic rotation has been used to orient the Wilson line in the temporal direction to reduce finite-volume effects at large z. The curve shows the perturbative result based on the analytic three-loop calculation [66,67], the analytic partial four-loop calculations [68][69][70][71][72][73][74], and the numerical full four-loop calculation [72]; its error band indicates the size of the O(α 4 s ) contribution.
where m is linearly divergent, Z φ is logarithmically divergent, and r mix is finite and associated with chiral symmetry breaking on the lattice. For our choices of Γ, the anticommutator vanishes and the expression simplifies to We follow the approach in Refs. [45,47] to determine m and Z 2 φ (1 − r 2 mix ), using the RI-xMOM intermediate scheme and converting to MS. Calculations are performed using the most chiral N f = 4 twisted mass ensembles from Ref. [59], averaging over pairs of ensembles with opposite PCAC masses rather than directly working at maximal twist. In addition to the operator with stout-smeared links used for the bare nucleon matrix elements, we also employ unsmeared links, which are expected to have reduced discretization effects, in some intermediate steps. After fixing to Landau gauge, we compute the position-space ζ propagator the momentum-space quark propagator, where χ is a quark field in the twisted basis; and the mixed-space Green's function for φ, These renormalize as To fix m, we evaluate the effective energy of the ζ propagator, which is renormalized by adding m. We use the nearest-neighbor lattice derivative. The relative matching among the three lattice spacings is done at z ≈ 0.61 fm. The absolute value of m is determined using unsmeared links on the finest lattice spacing, which is expected to produce the smallest discretization effects, and matching to the determined using unsmeared links in the RI-xMOM scheme, matched to MS, and evolved to scale 2 GeV. A hypercubic rotation has been used to orient the Wilson line and the quark momentum in the temporal direction. The multiple points at the same µ 2 have different values of the RI-xMOM scheme parameter y. The lines with error bands give the extrapolation to a 2 µ 2 = 0 and the darker part of each error band indicates the fit range. (1)(57) 0.986(25) 1.95 −0.373 (1)(50) 0.908(20) 2.10 −0.305 (1)(37) 0.907(11) TABLE IV. Renormalization parameters from the auxiliary-field approach, determined via the RI-xMOM intermediate scheme.
The second uncertainty for the auxiliary field mass comes from the absolute matching onto perturbation theory and is fully correlated across the three ensembles.
The other renormalization factors are determined using conditions designed to eliminate dependence on m: 1 12 These are evaluated at the scale µ 2 = p 2 , choosing p = p zẑ . This defines a family of renormalization schemes that depend on the dimensionless quantity y ≡ p z z. From the above, we extract the relevant overall renormalization factor, Z 2 φ (1 − r 2 mix ), at fixed kinematics. We then convert Z φ to the MS scheme using the one-loop expression from Refs. [45,47] and evolve to the scale 2 GeV using the two-loop anomalous dimension of the static-light current [75,76]. Fig. 7. As this is done at relatively high scales where the perturbative matching and evolution are applicable, we do this using unsmeared gauge links. Except at low µ 2 , the statistical uncertainty is negligible compared with systematics. At each µ 2 , we estimate the latter such that the spread of results for different scheme parameters y is covered. For each lattice spacing, we extrapolate a 2 µ 2 to zero assuming a linear dependence; the systematic uncertainty is propagated assuming a 50% correlation between every pair of points. Following the approach used in Ref. [45], we match between unsmeared and smeared links in the infrared regime at large z and small p 2 .
Final parameters for operators with stout-smeared links are given in Table IV. The large uncertainty for the mass parameter is caused by the absolute matching onto perturbation theory. At each z, this absolute matching produces an overall factor applied to h Γ (z) at all three lattice spacings. Therefore it can be ignored when studying the approach to the continuum limit. However, this uncertainty must be included when comparing continuum-limit results against other renormalization approaches.
An additional perturbative conversion [12] yields results in the MMS scheme; this cancels a log(z 2 ) divergence in the MS-renormalized matrix element at short distance. However, this conversion has only been computed at oneloop order, meaning that the cancellation may be inexact and some part of the divergence may still remain. The renormalized nucleon matrix elements for the three lattice spacings are shown in Fig. 10.

C. Ratio with zero-momentum matrix element
The simplest way to cancel ultraviolet divergences is to compute matrix elements of the same operator in different hadronic states and then take their ratio. Here we choose to take the ratio of matrix elements in a nucleon at nonzero momentum (i.e. those used throughout this paper) with the same in a nucleon at rest, As the signal-to-noise problem is much milder in a nucleon at rest, this requires a relatively inexpensive additional calculation: see Table V.
This ratio is similar to the reduced Ioffe-time distribution used in the pseudo-PDF approach for parton distributions [77]. Although it is a different observable than the MMS-renormalized matrix elements used for quasi-PDFs, it provides the opportunity to study the approach to the continuum limit in a clean, controlled setting. As such, this section can be seen as a preview of the continuum extrapolations of the renormalized matrix element h Γ in Section V A.
We consider variations of the continuum extrapolation in two different ways. First: precisely which points should be used to obtain R Γ (p z , z) at zero lattice spacing? One option is to ignore small differences in p z among the three ensembles, interpolate the lattice data to a common value of z in physical units, and then perform the extrapolation. Alternatively, noting that parameter x of quasi-PDFs is Fourier-conjugate to the product zp z , we can choose to interpolate to a common value of zp z before extrapolating; this could be more reliable because it accounts at leading order for the small differences in p z .
Second: what fit form should be used? As we have three lattice spacings, we restrict ourselves to two-parameter fits. At z = 0, the operator O is local, namely a vector or axial current; since we work at maximal twist, this calculation benefits from automatic O(a) improvement [78,79] and we extrapolate using an affine function of a 2 . For z = 0 (i.e. for the bulk of our data), O is nonlocal and there can be O(a) contributions that are not eliminated by automatic improvement [47]. In practice, it is not clear whether we are in the regime where O(a) contributions dominate; therefore, we extrapolate using both affine functions of a and of a 2 .
The ratio data and their extrapolations are shown in Figs. 8 and 9. When plotted versus z in physical units, clear discrepancies between the three ensembles are visible and for most of the parameter space, the lattice data from the coarsest ensemble are more than one standard deviation away from the extrapolations. These discrepancies are reduced when plotting the data versus zp z , although they remain significant for the unpolarized data at large z.
From this study, it appears that performing the extrapolation at fixed values of zp z is a better approach. For the unpolarized matrix elements with zp z < 5 and for the helicity matrix elements, lattice artifacts have a modest effect and are well under control, with the O(a) and O(a 2 ) extrapolations in good agreement. For the unpolarized matrix elements with zp z > 5, there is a stronger dependence on a and worse agreement between the two extrapolations; this suggests that at longer distances the lattice artifacts are less well controlled.

A. Renormalized matrix elements
Based on our study of the ratios of matrix elements in the previous section, we choose to linearly interpolate our MMS-scheme renormalized matrix elements to common values of zp z and then perform continuum extrapolations at each interpolated point. We again extrapolate in two ways, assuming lattice artifacts are either linear or quadratic in the lattice spacing.   Matrix elements renormalized using the whole-operator approach are shown in Fig. 11, along with their continuum extrapolations. Qualitatively, the picture is similar to the auxiliary-field renormalization approach, except that at small zp z , the real part of the matrix elements from the three lattice spacings are in better agreement, producing a milder effect from the continuum extrapolation and a better agreement between the two extrapolations. The latter is especially true for the unpolarized matrix element. Details of these continuum extrapolations for selected values of zp z are shown in Figs. 12 and 13. Clearly, our lever arm in a is limited, which makes it difficult to detect a preference for either of the two fits; this also produces a large uncertainty for the O(a) extrapolations.
Results from the two renormalization approaches are compared in Fig. 14. The whole operator approach tends to produce a smaller central value and a smaller uncertainty than the auxiliary field method. For the imaginary part of the matrix elements, the O(a 2 ) auxiliary-field extrapolation is in significant disagreement with both of the whole-operator extrapolations. In contrast, the O(a) result using the auxiliary field method is largely compatible with both whole-operator extrapolations, for low to medium values of zp z . This suggests that there may be significant O(a) lattice artifacts in the determination of the auxiliary field renormalization parameters and that it is necessary to account for them when taking the continuum limit.
Since renormalization in the auxiliary field approach is determined by just two parameters, one might ask whether there exist parameters that produce results compatible with the whole operator method. Figure 15 shows the effect of reducing the magnitude of the auxiliary-field mass renormalization parameter by δm = 0.4 GeV. Although this adjustment is hard to justify from the analysis in Section IV B, in Ref. [58] it was shown that its effect on quasi-PDFs  is suppressed by the factor δm/p z at large momentum. This change produces good agreement for the imaginary part of the matrix elements. However, some discrepancies remain for the real part, particularly in the unpolarized case at small zp z , where the slope of the auxiliary-field result is considerably steeper than the whole-operator data.
In the rest of this paper where we examine the effect on parton distributions, we will focus on the more precise data renormalized using the whole-operator method. However, we will continue to compare O(a) and O(a 2 ) extrapolations since they are not in complete agreement and we have no a priori reason to prefer one over the other.

Comparison with phenomenology
Before transforming the position-space matrix elements to obtain PDFs and comparing directly with phenomenology, we perform the reverse exercise. Starting with phenomenological parton distributions determined by NNPDF [80,81], we invert the matching and the Fourier transform to determine the position-space matrix elements that yield those PDFs, up to higher-order corrections in the matching. Figure 16 compares this with the continuum-extrapolated  [80,81] (dark gray). Note that in the lattice calculation, the pion mass is much larger than in nature, so that exact agreement should not be expected.
lattice matrix elements. Full agreement cannot be expected, since the lattice calculation was done at a heavy pion mass and other systematics such as the dependence on p z and finite-volume effects have not been included in this study.
The real part of the unpolarized matrix elements show reasonable agreement for zp z < 5; in the same range, the helicity matrix elements from the lattice lie below those from phenomenology. The helicity case can be partly understood by recalling that at heavy pion masses, the nucleon axial charge (i.e. the helicity matrix element at z = 0) lies below its physical value. At short distances, the imaginary parts of the lattice data have larger (more negative) slopes than phenomenology; the O(a) extrapolations are consistent with the latter at the 1σ level whereas the O(a 2 ) extrapolations are not. At nonzero lattice spacing, the slope is even larger and in worse agreement with NNPDF, so that the continuum extrapolation produces results that lie closer to phenomenology.
At larger values of zp z , there is a qualitative difference: the phenomenological curves tend steadily toward zero, whereas the lattice data do not. This is especially true for the unpolarized lattice matrix elements, of which both the real and imaginary parts are positive and increasing at large distances. At the coarsest lattice spacing, the lattice data lie well below zero (see Fig. 11), so it appears that the continuum extrapolation may be an overcorrection. Another way to characterize the imaginary part is via the position of the minimum of the curve: in the lattice data, it lies at a shorter distance than in phenomenology. This is consistent with the general expectation that correlation functions are shorter ranged at heavier pion masses.

B. Parton distributions
In this section, we present the main results of this paper, namely the effect of the continuum extrapolation on PDFs. However, we first discuss another source of systematic uncertainty: how to perform the Fourier transform in the definition of the quasi-PDF using a finite set of position-space data. We illustrate this using data on the finest ensemble, D45. Next, we perform the continuum extrapolation at fixed x, using the PDFs determined on each ensemble, and compare the result with the PDF determined from the continuum-limit matrix elements obtained in the previous section. Finally, we compare our continuum-limit PDFs with phenomenology.

Reconstruction techniques
As given in Eq. (3), the quasi-PDFq(x) is obtained from a Fourier transform (FT) of the renormalized matrix elements h(z). In practice, we obtain h(z) at intervals of the lattice spacing 2 , i.e. z/a ∈ Z. It is also necessary to truncate the FT at |z| ≤ z max , both because of the finite lattice size, which imposes z max < ∼ L/2, and because of growing statistical uncertainty at large |z|. Together, these have the effect of replacing the continuous FT by a truncated discrete FT (DFT): The discrete sampling makes the result formally periodic, so that it must be cut off at |x| ≤ π/(ap z ), which is at least 4 in our setup. The truncation introduces an additional systematic uncertainty [12], as shown using ensemble D45 in Fig. 17 for quasi-PDFs and PDFs. The latter are obtained by applying the matching procedure and nucleon mass corrections [17]. For the quasi-PDF, the effect of truncation is that one obtains a convolution of the desired result:q so that any features narrower in x than (z max p z ) −1 are smeared out. This is clearly visible in Fig. 17, where smaller values of z max p z are associated with broader quasi-distributions. The effect is reduced after applying the matching to obtain PDFs: results with z max p z = 4.7 and 5.9 are very similar. However, for z max p z = 3.5, both the unpolarized and helicity PDFs have qualitatively quite different behaviour, with a higher value for x between roughly −0.7 and −0.1 and a larger slope for x less than −0.3 as well as a smaller slope at small positive x and a peak at larger x in the positive region.
Since the Fourier transform introduces a systematic uncertainty, we supplement the naïve truncated FT with more sophisticated reconstruction techniques [82,83]. In these approaches, obtaining the Fourier transform from a finite number of data points is seen as an ill-defined inverse problem. Its solution is not unique and one approach is to use explicit models for the shape of the (quasi-)PDF. By contrast, we choose to use two approaches that do not contain an explicit model: the Backus-Gilbert method, first applied for PDFs calculations in [82] and the Bayes-Gauss-Fourier Transform (BGFT) [83]. These two procedures address the reconstruction problem as follows.

Backus-Gilbert (BG):
The inverse problem is obtained by inverting Eq. (3) to write the real and imaginary parts of the unpolarized matrix element in terms of the quasi-PDF: where for x ≥ 0,q ± (x) =q(x) ±q(−x), and likewise for the helicity case 3 . The reconstruction is applied independently toq + andq − , so for brevity we describe the procedure applied toq + . We also omit the labels p z and µ. For each x, the solution is assumed to be a linear combination of the finite set of lattice data: where a + can be understood as an approximation to the inverse of the Fourier transform in Eq. (25). The accuracy of this approximation is governed by the function that approximates δ(x − x ). Specifically, the result is an integral over the quasi-PDF: The function a + is determined by the Backus-Gilbert procedure [84], which minimizes the width of ∆ + (x, x ). For more details, see Refs. [40,82].
Bayes-Gauss-Fourier Transform (BGFT): Rather than directly reconstructingq from the lattice data, this procedure reconstructs a continuous form of the position-space matrix elements for all values of z: For this, we apply a nonparametric regression technique, based on Bayesian inference, called Gaussian process regression (GPR) [85]. This allows us to incorporate into the prior distribution the asymptotic behavior of the matrix elements (expected to decay to zero), as well as their smoothness properties. The result is continuous, defined for all real z, and has a Fourier transform computable in closed form. Taking the FT of h GPR (z), we refer to the result asq BGFT (x). More details are given in Ref. [83].
In Fig. 18, we compare results from the truncated discrete Fourier transform, Eq. (23), and the BG and BGFT reconstruction methods described above, again using ensemble D45 as our reference data set. For a fair comparison, in all cases we use z max p z = 4.7. We begin by discussing the quasi-PDFs (upper two panels). The most striking difference is that the Backus-Gilbert result has a discontinuity at x = 0 that is not present in the other results. This is becauseq BG − (x) is not constrained to vanish at x = 0. Such a discontinuity could occur if h(z) has a slowly decaying tail ∼ 1/z. For x between −0.5 and 1.0, the DFT and BGFT results are similar, although the BGFT distribution is slightly narrower. For larger values of |x|, the DFT produces stronger oscillations, which are suppressed by the BGFT. The BG result is the outlier, being considerably smaller at small negative x and also having a smaller dip below zero.
We next discuss the physically relevant parton distributions, obtained after matching and nucleon mass corrections (lower two panels). For most values of x, the DFT and BGFT method produce very similar results, although for BGFT the the dip below zero in the antiquark region occurs at smaller negative x and the magnitude is smaller at x = −1 and +1. Again, the BG result is somewhat different: in the antiquark region at small negative x, the small positive bump is gone and the result is either consistent with zero (unpolarized) or slightly negative (helicity). This discrepancy at small x may be associated with a lack of data for the matrix element at large |z|; better data or a more rigorous understanding of the large-|z| behavior could help to improve this situation. In the quark region for x greater than about 0.5, the BG result has a much weaker downward trend than the other two methods. Given that the DFT produces a result not substantially different from BGFT, we exclude the DFT from further analyses presented in the next sections.

Continuum extrapolation
In what follows, we compare the distributions at finite lattice spacings with continuum extrapolations. In the reconstruction of the quasi-PDFs we use the lattice data with |zp z | ≤ z max p z = 4.7, at which point either the real part or the imaginary part of the continuum matrix element is compatible with zero, as shown in Fig. 11. Moreover, we estimate the systematic uncertainty from this choice of the cutoff by varying z max : Finally, we estimate the combined uncertainty as the quadrature sum of ε cutoff (x) and the statistical uncertainty.
One approach for obtaining continuum-limit PDFs is to take the PDF determined on each ensemble and then perform an O(a) or O(a 2 ) extrapolation of the data at each x. This is shown in Fig. 19, for both unpolarized and helicity PDFs determined using the BG and BGFT methods. In the quark region with x between roughly 0 and 0.7, the PDFs decrease monotonically with the lattice spacing; at larger x, the D45 data (with the finest lattice spacing) move relatively upward to lie between those of the other two ensembles. For all x > 0, the O(a 2 ) extrapolation lies below all of the individual lattice spacings and the O(a) extrapolation is even lower. Using the BGFT approach, both of the extrapolations are consistent with the expected value of zero at x = 1, whereas for BG, this is true only of the O(a) extrapolation. In the antiquark region, the extrapolated results lie above the PDFs determined at finite lattice  Fig. 19) and at fixed zpz (pink, based on the continuum-limit data in Fig. 11). The distributions has been obtained using the BG (top panels) and BGFT (bottom panels) reconstruction techniques.
spacing, except for the BGFT unpolarized distribution near x = −1. This produces a more prominent positive region at small negative x, particularly in the unpolarized case. At larger negative x, the extrapolations are generally closer to zero.
Another approach is to obtain PDFs from the continuum limit of h(z) as determined in Section V A by extrapolating data at fixed zp z . By changing the order in which the continuum limit and the combination of the Fourier transform and PDF matching are performed, we obtain results affected by different systematic effects. The comparison of the  Fig. 20. They are consistent within uncertainties, except near x = 1, where the fixed-x extrapolation is in all cases lower than the fixed-zp z extrapolation and only the former is consistent with zero at x = 1.

O(a) extrapolations from both approaches is shown in
For comparing with phenomenology in the next section, we take the fixed-x extrapolation as our central value and add an additional systematic uncertainty in quadrature, namely half the difference with the fixed-zp z extrapolation.

Comparison with phenomenology
In Fig. 21, we compare the distributions obtained using O(a) and O(a 2 ) extrapolations with those obtained from phenomenology by NNPDF [80,81]. This comparison is intended to be qualitative, since our calculation was not done at the physical pion mass and does not include a study of other sources of systematic uncertainty such as finite-volume effects or the dependence on p z .
In the antiquark region (x < 0), the NNPDF result is slightly positive for x > −0.25, particularly in the unpolarized case. Focusing on the latter case, both of the extrapolations using both BG and BGFT methods reproduce this feature, although the O(a) extrapolation (which has a larger uncertainty) prefers a wider and larger positive region. This agreement with NNPDF is only present after the continuum extrapolation and does not appear in the analyses of any of the individual ensembles. For larger negative x, the NNPDF distributions are close to zero. However, the BGFT result is below zero, particularly when using an O(a 2 ) extrapolation.
In the quark region (x > 0), the distributions obtained from our data tend to have smaller peaks at larger x than phenomenology and fall off more slowly at large x. All of the analyses are consistent with zero at x = 1, except for the O(a 2 )-extrapolated BG data. For small x, the lattice unpolarized distributions are consistent with phenomenology, whereas the lattice helicity distributions have smaller slopes. In the unpolarized case, the agreement holds for a wider range of x when using the BGFT approach, and this approach also produces less disagreement in the helicity distribution.

VI. CONCLUSIONS
In this work we performed a lattice QCD calculation of isovector parton distributions via the quasi-PDF approach, using three twisted mass ensembles with different lattice spacings. This enabled a study of discretization effects, which can first appear at linear order in the lattice spacing, and the approach to the continuum limit.
Although our data are unable to clearly distinguish O(a) from O(a 2 ) contributions, we nevertheless observed significant discretization effects, both in the position-space matrix elements and in the final parton distributions. In the antiquark region, taking the continuum limit produces a reasonable agreement with phenomenology. Previous calculations, such as the one at the physical pion mass in Ref. [12], have failed to reproduce the phenomenological behaviour at small negative x; our work suggests that discretization effects contribute significantly to this discrepancy. At larger negative x, the agreement is better when using the Backus-Gilbert method, although the uncertainty is also larger. In the quark region, the continuum extrapolation also has a significant effect, although large disagreements with phenomenology remain. The latter is unsurprising, as we have not controlled other sources of systematic uncertainty.
Going beyond the naïve truncated discrete Fourier transform, we have compared two reconstruction techniques for obtaining quasi-PDFs from a finite set of lattice data. We found that the Bayes-Gauss-Fourier-Transform method produces a somewhat better agreement with phenomenology in the quark region and worse agreement in the antiquark region, although for the latter the Backus-Gilbert method has a larger uncertainty. Given the uncontrolled systematic effects, these observations should be treated with caution.
We have also compared two different approaches for nonperturbative renormalization of the nonlocal operator O Γ (z). The auxiliary-field approach tends to produce significantly larger renormalized matrix elements than the whole-operator approach, particularly at large z. In this work we chose to study PDFs using the latter because its results are more precise, but it will be important to continue studying different renormalization approaches to understand their different systematics and whether they all produce the same continuum limit.
While we have demonstrated the importance of discretization effects, more work will be needed to understand the relative importance of O(a) and O(a 2 ) effects in typical calculations. This could be done by performing calculations using a wider range of lattice spacings or by applying Symanzik improvement to remove O(a) effects [47].