Flavor non-singlet parton distribution functions from lattice QCD at physical quark masses via the pseudo-distribution approach
Manjunath Bhat, Krzysztof Cichy, Martha Constantinou, Aurora Scapellato
PParton distribution functionsfrom lattice QCD at physical quark massesvia the pseudo-distribution approach
Manjunath Bhat, Krzysztof Cichy, Martha Constantinou, and Aurora Scapellato Faculty of Physics, Adam Mickiewicz University,ul. Uniwersytetu Pozna´nskiego 2, 61-614 Pozna´n, Poland Temple University, 1925 N. 12th Street, Philadelphia, PA 19122-1801, USA (Dated: May 6, 2020)One of the great challenges of QCD is to determine the partonic structure of the nucleon fromfirst principles. In this work, we provide such a determination of the unpolarized parton distribu-tion function (PDF), utilizing the non-perturbative formulation of QCD on the lattice. We applyRadyushkin’s pseudo-distribution approach to lattice results obtained using simulations with thelight quark mass fixed to its physical value; this is the first ever attempt for this approach directlyat the physical point. The extracted coordinate-space matrix elements are used to find the relevantphysical Ioffe time distributions from a matching procedure. The full Bjorken- x dependence of PDFsis resolved using several reconstruction methods to tackle the ill-conditioned inverse problem en-countered when using discrete lattice data. Another novelty of this calculation is the considerationof the combination with antiquarks q v + 2¯ q . The latter, together with the non-singlet valence quarkPDF q v , provides information on the full distribution. Good agreement is found with PDFs fromglobal fits already within statistical uncertainties and it is further improved by quantifying severalsystematic effects. The results presented here are the first ever ab initio determinations of PDFsfully consistent with global fits in the whole x -range. Thus, they pave the way to investigating awider class of partonic distributions, such as e.g. singlet PDFs and generalized parton distributions.Therefore, essential and yet missing first-principle insights can be achieved, complementing the richexperimental programs dedicated to the structure of the nucleon. I. INTRODUCTION
Despite the fact that the nucleon is the main buildingblock of visible matter and is responsible for almost allof the mass of the visible Universe, it is only now thatseveral aspects of its internal structure are beginning tobe thoroughly explored. The wealth of data from present-day experiments, e.g. from the Large Hadron Collider andthe Jefferson Laboratory 6 and 12 GeV programs, allowsus to unravel many details that so far eluded any insight.Moreover, the planned and recently approved Electron-Ion Collider at the Brookhaven National Laboratory, thatwill start operation in around ten years, is oriented at an-swering important questions about the nucleon structure,such as of the origin of the proton mass, the spin distribu-tion and the role of gluons [1]. Along with the experimen-tal efforts, there is constant progress in the theoreticalunderstanding, based on exhaustive analyses of empiri-cal data and various approaches to describe the physics ofthe strong dynamics of partons (quarks and gluons), gov-erned by the theory of quantum chromodynamics (QCD).One of the most important tools is perturbation the-ory, which is, however, limited to high energy scales,at which the perturbative expansion can evince conver-gence. Meanwhile, it is clear that all energy scales con-tribute to the dynamics of the nucleon and hence thedescription of the non-perturbative aspects is of essentialimportance. This can take the form of phenomenologicalmodels, which have led to important successes. However,a truly ab initio knowledge is, in principle, still possibleto extract directly from the QCD Lagrangian. The most successful non-perturbative approach to QCD is to for-mulate it on a discrete spacetime grid, the lattice. Thisleads to a regularization of the QCD path integral andresults in multidimensional integrals that can be evalu-ated numerically, usually with Monte Carlo simulations.Such computations, however, are necessarily done in Eu-clidean spacetime, which poses a fundamental problemfor partonic physics, most naturally formulated in termsof light-cone correlations. The latter cannot be accessedin Euclidean lattice QCD (LQCD) and thus, the informa-tion on the nucleon structure from LQCD has been lim-ited for many years. Most of the insights were obtainedfrom lattice calculations of moments of parton distribu-tion functions (PDFs) and generalized parton distribu-tions (GPDs). In principle, the full distributions can bereconstructed from a sufficiently large number of theirmoments, as obtained from LQCD. However, in practice,the computations are limited to only the lowest 2-3 mo-ments. For higher moments, the breaking of rotationalsymmetry implied by the lattice leads to unavoidablepower-divergent mixings with lower-dimensional opera-tors and what is more, the signal-to-noise ratio for suchhigher moments is quickly decaying. For recent calcula-tions of the low moments by the Extended Twisted MassCollaboration (ETMC), we refer the interested reader toRefs. [2–7].Extending the calculations of moments to the full x -dependence of the partonic distributions has been a sub-ject of intense studies over the years. Even though theearly proposals date back to the previous century, the re-cent revival of this topic, one that brought considerable a r X i v : . [ h e p - l a t ] M a y progress to this field, came with the proposal of Ji [8] tocalculate so-called quasi-distributions. They are relatedto the desired light-cone distributions, but computableon a Euclidean lattice. The underlying idea is to replacethe light-cone correlations in the definition of a distri-bution by spatial ones between boosted nucleon states.One can then exploit the fact that the obtained quasi-distributions share the infrared physics with their light-cone counterparts. As a consequence, their difference isin the ultraviolet and can be matched using perturbationtheory, utilizing the so-called Large Momentum EffectiveTheory (LaMET) [9]. Ji’s proposal sparked a huge tideof theoretical and numerical efforts to understand crucialaspects of this approach, such as renormalizability andappropriate renormalization prescriptions, matching, nu-cleon mass corrections, higher-twist effects, finite volumeeffects, as well as to extract the distributions in varioussetups, see e.g. Refs. [10–54].Quasi-distributions can be thought of as generalizationof the light-cone ones to the finite-momentum frame. Asshown by Radyushkin in a series of papers [36, 55–61], thesame matrix elements that specify a quasi-distributioncan also be used to define another generalization of itslight-cone counterpart, the so-called pseudo-distribution.These matrix elements can be viewed as functions of twoLorentz invariants. The first one is the spacetime in-terval z , where z µ can be taken as (0 , , , z ) and de-scribes the separation between the quarks in the insertedoperator (containing the Wilson line to guarantee gaugeinvariance). The other Lorentz invariant is the product ν ≡ − p · z (with p µ being the nucleon boost) and the vari-able ν is called the Ioffe time [62]. The matrix elementswritten as the function of z and ν are called Ioffe-timedistributions (ITDs).The underlying difference between quasi- and pseudo-distributions is that the former are defined as the Fouriertransform of the matrix elements in z , while to obtainthe latter, one takes the transform in the Ioffe time.This has far-reaching consequences and makes the ap-proaches inequivalent, even though they can be com-puted from the same matrix elements. In particular,pseudo-distributions have the canonical support in theBjorken- x fraction, − ≤ x ≤
1, as opposed to quasiones that can be non-zero also for x outside of this range.The matching of pseudo-distributions to the light-coneframe [31, 36, 59, 63] is peformed at the level of ITDs.As we shall see below, the matching is numerically asmaller effect than for quasi-PDFs. In particular, it de-pends less significantly on the region of x . It is also animportant difference that pseudo-distributions can fullyutilize lattice data at all nucleon boosts, i.e. they con-tain physical information also from low momenta, effec-tively at small Ioffe times. Thus, the matching of ITDs isnot based on LaMET, but still on a factorization of thefinite-momentum distribution into a light-cone one anda perturbatively computable matching coefficient. Thepseudo-distribution approach to extract the valence un-polarized PDF of the nucleon was explored first in the quenched setup [64] and recently it was extended to in-clude the effects of dynamical quarks at non-physicalmasses [65] and to investigate the PDFs of the pion [66].Even more recently, it was also applied in a setup withlight quark mass close to, but not yet directly at its phys-ical value, with the corresponding pion mass of 170 MeV[67]. In addition, moments of ITDs were computed [37]and the general issue of reconstructing distributions fromITDs under incomplete Fourier transforms was analyzedin related studies by the same group [44].Apart from approaches based on ITDs, we also men-tion other proposed methods that can lead to determina-tions of the full x -dependence of partonic distributions.Instead of using a Wilson line to ensure gauge invari-ance, they can employ auxiliary propagators, of fictitiousscalar [68], heavy [69] or light quarks [70–72]. Approachesbased on the hadronic tensor also exist [73, 74]. Simi-larly to quasi- and pseudo-distributions, these methodsare also intensely investigated by various groups, see e.g.Refs. [75–78]. An extensive review of these efforts, withemphasis on the most investigated approach of quasi-distributions, can be found in Ref. [79]. Recently, a re-view devoted to the principles and various applicationsof LaMET also appeared [80].In this paper, we apply the pseudo-distribution ap-proach for the first time to lattice data obtained withlight quark mass fixed to its physical value. Anothernovel aspect is the extraction of not only the valencedistribution, denoted here by q v ( x ), but also of the com-bination with antiquarks q v ( x ) + 2¯ q ( x ). In this way, weare able to present results for the full distribution, q ( x ) = q v ( x ) + ¯ q ( x ) and for the sea quark PDF, q s ( x ) = ¯ q ( x ).We explore several systematic effects inherited in latticecomputations and we address the issue of reconstructingthe PDFs from ITDs, subject to an ill-defined inverseproblem. We show that this leads ultimately to full con-sistency between all considered distributions and the cor-responding ones from global fits, in the whole range ofBjorken- x . Even though at the present stage full quan-tification of all systematics is not yet possible, the resultsobtained in this work are unambiguously optimistic anddemonstrate the potential of the techniques for an issuethat was for many years thought to be too difficult forthe lattice.The outline of the remainder of the paper is the fol-lowing. In Sec. II, we discuss theoretical principles ofpseudo-PDFs and practical aspects of their computa-tions. Then, we present the lattice details of the cal-culation in Sec. III. In Sec. IV, we show our results forthe pseudo-distributions and their matching to light-coneITDs and we compare different methods for the recon-struction of PDFs. We also discuss systematic effectsfrom the choice of the Ioffe time range, the value of thestrong coupling constant and from other sources and inSec. IV D, we show our final PDFs. Finally, we concludeand discuss future prospects in Sec. V. II. THEORETICAL SETUP AND ANALYSISTECHNIQUES
We start by summarizing the relevant steps in theprocedure leading from lattice-extracted matrix elements(pseudo-ITDs) to the light-cone distributions, ITDs andPDFs after a suitable reconstruction. We refer theReader to the review of Ref. [60] for an extensive dis-cussion on the theoretical principles and properties ofpseudo-distributions.The underlying matrix elements computed on the lat-tice are defined (in Euclidean spacetime) as M ( ν, z ) = (cid:104) P | ψ (0 , z ) γ W ( z, ψ (0 , | P (cid:105) , (1)where | P (cid:105) is a boosted nucleon state with four-momentum P µ =( P , , , P ) and W ( z,
0) is a straightWilson line. The Wilson line is chosen along the directionof the boost and has length z (i.e. we take z µ = (0 , , , z )and z will henceforth refer to the length of z µ ). With thischoice of kinematics, the Ioffe time ν = P z . The Diracstructure γ leads to a faster convergence to the light-cone ITD, as compared to another plausible choice of the γ structure [55]. It was also found that γ avoids a finitemixing with the twist-3 scalar operator due to the break-ing of chiral symmetry by some fermionic discretizations[23].The matrix element defined by Eq. (1) exhibits twokinds of divergences: standard logarithmic one and apower divergence related to the Wilson line. However,it has been shown to be multiplicatively renormalizableto all orders in perturbation theory [21, 22]. In Ref. [64],it was suggested that the divergences can be canceled byforming a double ratio with zero-momentum and local( z = 0) matrix elements: M ( ν, z ) = M ( ν, z ) / M ( ν, M (0 , z ) / M (0 , . (2)We follow this renormalization procedure and we refer to M ( ν, z ) as reduced matrix elements or pseudo-ITDs. Itis important to emphasize that the double ratio not onlyremoves the divergences, but also part of the higher-twistcontamination, which is generically of O ( z Λ ) [64].The double ratio defines a renormalization scheme wherethe renormalization scale is proportional to the inverselength of the Wilson line.To get from pseudo-ITDs to light-cone ITDs and fi-nally to light-cone PDFs, a matching procedure is re-quired, similiarly to the quasi-PDFs case. The reducedmatrix elements, defined at different scales 1 /z , need tobe evolved to a common scale, 1 /z (cid:48) , and it is desirableto also convert them to the renormalization scheme com-monly used for PDFs, the MS scheme, where its renor-malization scale will be denoted by µ . The full matchingequation, to one-loop order in perturbation theory, reads [31, 36, 59, 63]: M ( ν, z ) = Q ( ν, µ ) + α s C F π (cid:90) du (3) × (cid:20) ln (cid:18) z µ e γ E +1 (cid:19) B ( u ) + L ( u ) (cid:21) Q ( uν, µ ) , where Q ( ν, µ ) is the MS-scheme light-cone ITD and thefunctions convoluted with Q are B ( u ) = (cid:20) u − u (cid:21) + , (4) L ( u ) = (cid:20) − u )1 − u − − u ) (cid:21) + (5)with the plus prescription defined as (cid:90) [ f ( u )] + Q ( uν ) = (cid:90) f ( u ) ( Q ( uν ) − Q ( ν )) . (6)The matching equation consists of two parts. The partcontaining the kernel B ( u ) evolves the pseudo-ITDs to acommon scale µ and the part with L ( u ) converts to theMS scheme.We invert the matching equation and to look sepa-rately into the effect of evolution and scheme conversion,we introduce intermediate evolved ITDs, M (cid:48) ( ν, z , µ ).Thus, M (cid:48) ( ν, z , µ ) = M ( ν, z ) − α s C F π (cid:90) du (7) × ln (cid:18) z µ e γ E +1 (cid:19) B ( u ) M ( uν, z ) . The evolved ITD has three arguments, the Ioffe time ν ,the common scale µ and the initial scale z . In principle,its value should be independent of the initial scale andwe will test this conjecture up to our statistical precision.The scheme conversion then follows: Q ( ν, z , µ ) = M (cid:48) ( ν, z , µ ) − α s C F π (cid:90) duL ( u ) M ( uν, z ) , (8)where again we will test the independence on the initialscale. For the reconstruction of the final PDF, discussedbelow, we will average the matched ITDs Q ( ν, z , µ ) forcases where a given Ioffe time is achieved by differentcombinations of ( P , z ) and denote such an average by Q ( ν, µ ).The ITDs, Q ( ν, µ ), are related to PDFs, q ( x, µ ), bya Fourier transform in Ioffe time: Q ( ν, µ ) = (cid:90) − dx e iνx q ( x, µ ) , (9)where the antiquark distribution for positive x is ¯ q ( x ) = − q ( − x ). Decomposing into real and imaginary parts andusing this property, one obtainsRe Q ( ν, µ ) = (cid:90) dx cos( νx ) (cid:0) q ( x, µ ) − ¯ q ( x, µ ) (cid:1) = (cid:90) dx cos( νx ) q v ( x, µ ) , (10)Im Q ( ν, µ ) = (cid:90) dx sin( νx ) (cid:0) q ( x, µ ) + ¯ q ( x, µ ) (cid:1) = (cid:90) dx sin( νx ) q v s ( x, µ ) , (11)which relate the valence distribution, q v = q − ¯ q , to thereal part of the ITDs and the other non-singlet distribu-tion involving two flavors, q v s ≡ q v + 2¯ q = q + ¯ q , to theimaginary part of the ITDs.All of the equations (9)-(11) involve a known left-handside (reduced ITDs computed on the lattice and subjectto the matching procedure) and integrals of a partonicdistribution to be determined. As discussed in detail inRef. [37], such determination poses an inverse problem,related to the fact that inverse equations are ill-defined.Namely, they involve an integral over continuous Ioffetime up to infinity, while on the lattice, one is necessarilyrestricted to a finite number of determinations of Q ( ν )that cover only a finite range of Ioffe time, from 0 up tosome ν max . To reconstruct the distributions, we will fol-low three ways. The inverse problem stems from havingincomplete information and hence, solving it is not possi-ble without additional assumptions. There is an infinitenumber of possible assumptions to provide the missinginformation and thus, it is important to use as mild onesas possible, in order not to bias the reconstruction pro-cedure. We will perform naive Fourier transforms andwe will use two additional ways of handling the inverseproblem. First, we will apply the Backus-Gilbert (BG)method [81], originally proposed to be used in PDF re-construction in Ref. [37]. Second, we will perform recon-struction by fitting the matrix elements using a fittingansatz for the light-cone PDF, as suggested in Ref. [65].The BG method minimizes the variance of the solutionto the inverse problem, i.e. it maximizes its stability withrespect to variation of the data within their errors. Thisvariance minimization condition is a model-independentassumption that provides a unique distribution given aset of input ITDs. For each value of Bjorken- x , the mini-mization condition defines a d -dimensional vector a K ( x )(where d is the number of available evaluations of the in-put ITD), which is an approximate inverse of the kernelfunction K ( x ) (the cosine or the sine function, respec-tively for Eqs. (10)-(11)), i.e.∆( x − x (cid:48) ) = (cid:88) ν a K ( x ) ν K ( x (cid:48) ) ν , (12)where K ( x (cid:48) ) is taken as a d -dimensional vector with el-ements K ( x (cid:48) ) ν = cos( νx (cid:48) ) or K ( x (cid:48) ) ν = sin( νx (cid:48) ). In theideal case of d → ∞ evaluations of the ITD spanning the infinite range of Ioffe times, thus defined vector a K ( x )would lead to ∆( x − x (cid:48) ) being the Dirac delta function δ ( x − x (cid:48) ). In turn, for finite d , a K ( x ) leads to an ap-proximation to the Dirac delta function with minimizedwidth. We find the vectors a K ( x ) from width minimiza-tion conditions, spelled out explicitly e.g. in Ref. [37],which yield a K ( x ) = M − K ( x ) u K u TK M − K ( x ) u K , (13)where the elements of the d × d -dimensional matrix M K ( x ) are given by M K ( x ) νν (cid:48) = (cid:90) dx (cid:48) ( x − x (cid:48) ) K ( x (cid:48) ) ν K ( x (cid:48) ) ν (cid:48) + ρ δ νν (cid:48) (14)and of the d -dimensional vector u K by u Kν = (cid:90) dx (cid:48) K ( x (cid:48) ) ν . (15)The parameter ρ in M K ( x ) regularizes this matrix(Tikhonov regularization [82], see also Refs. [37, 83, 84]),making it invertible. The value of ρ should be relativelysmall in order not to bias the result and not to decreasethe resolution of the method. In our study, we find that ρ = 10 − is a good compromise, with smaller values in-troducing large oscillations in the final distributions dueto the presence of very small eigenvalues of M K ( x ).Having found the vectors a K ( x ) (for both kernel func-tions), the distributions q v or q v s are reconstructed as q v/v s ( x, µ ) = (cid:88) ν a K ( x ) ν Re / Im Q ( ν, µ ) . (16)We will also consider a version of the BG procedure withpreconditioning (BG precond ), which can be applied in acase where a realistic guess of the solution of the inverseproblem is available. In this variant, the kernel functionand the desired distribution are rescaled by a function p ( x ), ˜ K ( x ) = K ( x ) p ( x ) and ˜ q ( x ) = q ( x ) /p ( x ). Then,˜ q ( x ) embodies the deviation of the reconstructed distri-bution from p ( x ). The procedure of calculating the vec-tors a K ( x ) is unchanged, apart from using the rescaledkernel and taking into account the preconditioning func-tion in the final reconstruction equation: q v/v s ( x, µ ) = (cid:88) ν a K ( x ) ν p ( x ) Re / Im Q ( ν, µ ) . (17)The other reconstruction technique that we will use isto assume a functional form of a fitting ansatz for thelight-cone PDF. This is analogous to procedures used inphenomenological fits of PDFs from experimental data.We will adhere to the simplest reasonable functional formthat captures the expected limiting behaviors at low andlarge- x for both q v and q v s , in the range x ∈ (0 , q ( x ) = N x a (1 − x ) b , (18)where the exponents a, b are fitting parameters. N isfixed to 1 /B ( a +1 , b +1) for q v , where B ( x, y ) is the Eulerbeta function, related to the gamma function, B ( x, y ) =Γ( x )Γ( y ) / Γ( x + y ). This ensures the normalization of thevalence distribution to 1. In the case of q v + 2¯ q , N is leftas an additional fitting parameter.The fits are performed minimizing the χ function de-fined as χ = ν max (cid:88) ν =0 Q ( ν, µ ) − Q f ( ν, µ ) σ Q ( ν, µ ) , (19)where σ Q ( ν, µ ) is the statistical error of the light-coneITD Q ( ν, µ ). Q f ( ν, µ ) is given by the cosine/sineFourier transform of the fitting ansatz (18), respectivelyfor fits of the real/imaginary part of ITDs. The fittingfunction is continuous and thus, this Fourier transform isnot subject to any inverse problem. We will refer to thevalues of Q f as “fitted” ITDs and they are a continuousfunction of the Ioffe time. Note the fits depend on themaximum Ioffe time, ν max , and we will investigate differ-ent choices of this parameter and the sensitivity of thefinal PDF results to this choice. III. LATTICE SETUP
The underlying matrix elements are the same as theones used for the computation of quasi-PDFs. Thus,we use the matrix elements of Eq. (1), calculated inRefs. [32, 45], corresponding to the Dirac structure γ ofthe unpolarized PDF case. For quasi-PDFs, the Fouriertransform defining the distributions is performed at afixed nucleon boost P and the data of Refs. [32, 45] con-cern the cases of P = 6 π/L, π/L, π/L , correspondingto 0.83 GeV, 1.11 GeV and 1.38 GeV in physical units,respectively. In the pseudo-PDF approach, the Fouriertransform is taken in Ioffe time and hence can profit alsofrom data at low nucleon momenta. Moreover, the doubleratio that defines the reduced ITDs requires the knowl-edge of the zero-boost matrix elements. Thus, for thiswork, we computed also the cases of P = 0 , π/L, π/L ,i.e. 0, 0.28 GeV and 0.55 GeV.For details of the computational techniques, we refer tothe broad description in Ref. [45]. Here, we summarizethe main aspects. We use one ensemble of gauge fieldconfigurations with two degenerate flavors of maximallytwisted mass fermions [85, 86] with a clover improvement[87], generated by the Extended Twisted Mass Collabora-tion (ETMC) [88]. The gauge action is Iwasaki-improved[89]. The bare quark mass was tuned to approximatelyreproduce the physical value of the pion mass ( m π =130 . m N = 932(4)MeV) [90]. The lattice spacing is a = 0 . ×
96 sites, which corre-sponds to a physical lattice extent of L ≈ . O ( a ) improvement of physical observables. How- ever, the matrix elements of Eq. (1) are not in this cate-gory, apart from the local ( z = 0) case. Thus, in general,the PDFs calculated in this work have leading cut-off ef-fects linear in the lattice spacing. As shown in Ref. [52],maximal twist can remove some of the O ( a ) contribu-tions, but still an explicit specific improvement programis necessary to fully eliminate them. Likewise, an im-provement program is needed also for other lattice dis-cretizations, including ones preserving chiral symmetry[52]. P P [GeV] N confs N meas π/L π/L π/L π/L π/L N confs )and the number of measurements ( N meas ) for each value ofthe nucleon boost used in this work. In Tab. I, we summarize the statistics for the com-putation of matrix elements. While the statistics for P < .
83 GeV can be further increased, we keep it lowerthan the three highest momenta. This is desirable, aswe aimed at a similar statistical precision for all data.At these low momenta, the signal-to-noise ratio is veryfavorable and thus, similar precision could be achievedalready with the modest number of measurements re-ported in Tab. I. For all momenta, we use a source-sinkseparation, t s , of 12 lattice spacings (1.13 fm), which wasshown [45] to suppress excited states effects to below thestatistical precision of the data. Since the contaminationfrom excited states increases at larger momenta, the ad-ditional matrix element computations for this work canbe safely assumed to be free from such effects when using t s = 12 a , at our level of statistical precision. IV. RESULTSA. Bare and reduced matrix elements
We start by showing bare matrix elements as a func-tion of z/a for the 6 computed nucleon boosts, from 0to 1.38 GeV, see Fig. 1. We observe that as the hadronmomentum is increased, the distribution of the real partof the matrix elements becomes slightly narrower withrespect to the length of the Wilson line, i.e. they decayto zero at smaller values of z . At the same time, theimaginary part becomes more pronounced for larger mo-menta. The local matrix element ( z = 0) is real andcontains no divergence. With the employed definition ofthe vector current, it is only subject to a normalizationfactor computed in Ref. [91], Z V = 0 . M (0 ,
0) is compatible with 1for all nucleon boosts.
FIG. 1. Real (top) and imaginary (bottom) part of the barematrix elements ( M ( ν, z ) at fixed P ) for the unpolarizedPDF. Shown are all nucleon boosts: P = 0 (green circles), P = 2 π/L (blue rhombuses), P = 4 π/L (red 5-stars), P =6 π/L (yellow 6-stars), P = 8 π/L (purple asterisks), P =10 π/L (cyan squares). We form the reduced matrix elements according toEq. (2) and plot them against the Ioffe time in Fig. 2,separately for each nucleon boost and for lengths of theWilson line z/a ∈ [0 , z/a smaller than approx. 8, reduced ITDs obtained fromdifferent combinations of ( P , z ) that lead to the sameIoffe time P z are compatible with each other within ourstatistical uncertainties. This suggests that as long asthe difference in the scale 1 /z is not too large, the scale-dependence of reduced ITDs is rather small. The clear-est deviations are observed in the real part for the lowestboost, P = 2 π/L , for z/a >
8, where the scales 1 /z correspond to around 250 MeV and below. B. Evolved and matched ITDs
The 1 /z scale-dependence of the ITDs at a fixed Ioffetime is accounted for by the evolution equation (7). Wechoose to evolve all matrix elements to the scale corre-sponding to µ = 2 GeV, which will become our finalrenormalization scale after scheme conversion to the MSscheme. Inspection of the logarithm in Eq. (7) revealsthat the scale µ = 2 GeV corresponds to ITDs beingevolved to z = 2 e − γ E − / /µ ≈ . a , i.e. around 0.067fm at our lattice spacing (1 /z ≈ . µ = 2 GeV, -101 0 2 4 6 800.511.5 FIG. 2. Real (top) and imaginary (bottom) part of the re-duced elements ( M ( ν, z ) at fixed P ) for the unpolarizedPDF. Symbols are the same as used in Fig. 1. α s /π ≈ . α s /π = 0 .
1, the latter taken inthe quenched study of Ref. [64] and close to the value ofRef. [65], α s /π ≈ . P , z ) at the same ν up to z/a = 10 in the real part andfor all values of z/a in the imaginary part.The final step of the procedure to arrive at light-coneITDs is to perform the scheme conversion according toEq. (8). The outcome is shown in Fig. 4. The agreementbetween data at a given Ioffe time coming from differentmomenta and lengths of the Wilson line holds, similarlyas for evolved ITDs, up to around 9-10 lattice spacingsin the real part and for all values of z in the imaginarypart.The effects of evolution and scheme conversion aresummarized in Fig. 5, where we averaged ITDs corre-sponding to the same Ioffe time ν , but originating fromdifferent combinations of ( P , z ). To avoid contamina-tion from ITDs that are off from a universal curve, werestricted the average to cases with z/a ≤
8. The ef-fects of evolution and scheme conversion are oppositeto each other and, interestingly, approximately equal inmagnitude. Thus, the final matched ITDs turn out to becompatible, within our statistical precision, with originalreduced matrix elements. It implies that the one-loopmatching procedure at the level of ITDs is a small effect,which raises hope that higher-order matching effects are
FIG. 3. Real (top) and imaginary (bottom) part of theevolved matrix elements ( M (cid:48) ( ν, z , µ ) at fixed P ) forthe unpolarized PDF. The scale after evolution is 1 /z = µe γ E +1 / / ≈ . µ = 2GeV. Symbols are the same as used in Fig. 1. -0.500.511.5 0 2 4 6 800.511.5 FIG. 4. Real (top) and imaginary (bottom) part of thematched MS( µ = 2 GeV) matrix elements ( Q ( ν, z , µ ) atfixed P ) for the unpolarized PDF. Symbols are the same asused in Fig. 1. even smaller. Nevertheless, obviously, a two-loop compu-tation is still desired to check explicitly this statement.It is also worth to contrast the one-loop matching for thecase of pseudo-PDFs with the one for quasi-PDFs. Inthe latter case, as e.g. in Ref. [45] that used the same FIG. 5. Real (top) and imaginary (bottom) part of the re-duced ( M ( ν, z ); orange circles), evolved ( M (cid:48) ( ν, µ ); greensquares) and matched ( Q ( ν, µ ); purple stars) Ioffe time dis-tributions. The matrix elements corresponding to the sameIoffe time ν coming from different combinations of ( P , z ) wereaveraged, keeping only the ones at z/a ≤ bare matrix elements as the present study, the one-loopmatching effects are considerably larger. The matchingfor quasi-PDFs is performed in x -space and the differ-ence between a quasi-PDF and a matched PDF are above100% in many regions of x . Hence, it is plausible thatthe matching in the pseudo-distribution approach, at thelevel of ITDs (in ν -space), is more controlled, i.e. lesssubject to truncation effects. C. Light-cone PDFs
We now present results for the unpolarized PDFs ob-tained from matched ITDs discussed in the previous sub-section.Reconstruction of a PDF from ITDs requires, in prin-ciple, the knowledge of the full Ioffe time dependence ofthe ITDs, from ν = 0 to ν = ∞ . Obviously, with numer-ical calculations of ITDs on the lattice, the upper limit, ν max , is necessarily finite. It is desirable to take ν max as large as possible, ideally to observe that ITDs havedecayed to zero. However, as Fig. 4 suggests, this is diffi-cult with the currently attained nucleon boosts. The realpart of matched ITDs approaches zero at ν ≈ −
8, butthese Ioffe times are obtained with Wilson line lengthsof order 1 fm at the largest boost, which correspondsto very low energy scales at which the matching proce-dure is likely to fail. For the imaginary part of ITDs,non-zero values are observed even at ν = 8. It is, thus,clear that a robust extraction of the full x -dependence re- FIG. 6. Top row, left/right panels: real/imaginary part of matched ITDs (blue circles / red squares) and fitted ITDs (blue/redband). Middle and bottom row: light-cone distributions q v (middle left), q v s (middle right), q (bottom left), ¯ q (bottom right)from 3 reconstruction methods: naive Fourier transform (orange), Backus-Gilbert (green), fitting ansatz for the distribution(purple). Shown are also NNPDF phenomenological distributions (grey) [92]. The ITDs are fitted up to ν max ≈ . z max /a = 4), α s /π ≈ . µ = 2 GeV. quires achieving even larger Ioffe times, especially for theimaginary part that yields the distribution q v s . Sincethe length of the Wilson line is limited by the reliabil-ity of the matching procedure and higher-twist effectsof O ( z Λ ), larger Ioffe times need to be reached atlarger nucleon boosts. This is, however, difficult for thelattice, which is due to the decaying signal-to-noise ra-tio with increasing nucleon momentum, as discussed inRef. [45].Meanwhile, we are in position to investigate what thecurrently available range of Ioffe times implies for the x -dependence. The key parameter to decide in reconstruct-ing PDFs is the maximum Ioffe time, ν max . Below, weprovide reconstructed PDFs for various choices of ν max ,ranging from around 2.6 to 7.9. The former correspondsto taking ITDs obtained from matrix elements with in-sertions of the operator with Wilson line length up to z max /a = 4 lattice units (around 0.37 fm) and the lat-ter to 12 lattice units (1.12 fm). It is unclear, a priori,which value of ν max ensures reliable matching and goodcontrol over higher-twist effects. However, given the sta- tistical uncertainties of our data, we adopt a criterionthat safe values of ν max are those for which the matchedITDs obtained from different combinations ( P , z ) corre-sponding to the same ν are consistent with each other.This criterion leads to maximum z of order 8-9 latticeunits (approx. 0.8 fm), as we have discussed in the con-text of Fig. 4. In this way, the reached Ioffe times are oforder 5-6. At these values of ν , the real part of matchedITDs is already close to 0, thus giving good hope for thereconstruction of the valence distribution. In turn, theimaginary part of ITDs is still rather far away from zero,which is expected to bring significant uncertainties intothe low- x behavior of q v + 2¯ q .In Figs. 6, 7 and 8, we show the reconstructed PDFswith z max /a = 4, 8 and 12, respectively ( ν max ≈ .
6, 5 . . z max /a = 4 (Fig. 6), the data cover an in-sufficient range of Ioffe times and thus, the naive Fouriertransform, as well as the BG method are simply missingthe data and lead to unrealistically looking distributions. FIG. 7. Same as Fig. 6, but the range of Ioffe times extended to ν max ≈ . z max /a = 8). However, interestingly, the PDF fitting ansatz approachprovides significantly better distributions. The data atsmall Ioffe times are very precise and guide the fits, lead-ing to precise results at large x , in full agreement withthe phenomenological curves of NNPDF [92] for all dis-tributions. The fitted matrix elements, Q f ( ν, µ ), aredepicted as bands in the upper row of Fig. 6 and the fitsprovide good description of the data ( χ / dof ≈ . . ν manifest themselves in the increasing uncer-tainty of the reconstructed PDFs at low x , in particularfor the q v s distribution coming from the imaginary partof matched ITDs. The latter uncertainty propagates it-self also to the full distribution q = q v + ¯ q and to thesea distribution ¯ q , obtained from linear combinations of q v and q v s . We note that the full distribution q is con-sistent with NNPDF for all values of x , which, however,holds for x < ∼ . x , but this difference is compensated in q by the perfect agreement of q v s in this range.We want to see now how robust are the PDFs ob-tained with ν max ≈ . z max /a = 8,which leads to ν max ≈ .
2. This range of Ioffe times con-tains significantly more data points than the range upto ν max ≈ . ν enteringthe fits ( χ / dof ≈ . .
4) for the fit of the real (imag-inary) part) provide more constraints for fitting param-eters and their most significant effect is to decrease thePDFs uncertainty at x < ∼ . − .
5. This originates fromthe reduced uncertainty in the fitted matrix elements, i.e.the smaller widths of bands in the upper row of Fig. 7 inthe region ν ≈ . − .
2, now constrained by the actuallattice data and not simply guided by the low- ν behavior.However, the bands overlap for the cases of z max /a = 4and z max /a = 8 and thus, the resulting PDFs move onlywithin the uncertainties of the former case. This is en-couraging, since it suggests relatively little dependenceon ν max , with the effect of increasing the latter restrictedpredominantly to giving more precise access to lower x values for the distributions.Similar conclusions are drawn when further increasingthe range of Ioffe times, to ν max ≈ . z max /a = 12). For0 FIG. 8. Same as Fig. 6, but the range of Ioffe times extended to ν max ≈ . z max /a = 12). this case, shown in Fig. 8, one needs to keep in mind thatthe contamination of the large- z ITDs might be signifi-cant, as best seen in the deviation of the z/a > ∼
10 lowest-momentum points from an universal curve (Fig. 4). Theextended ν -range has a similar effect on the distribu-tions as the one when increasing ν max from 2.6 to 5.2– PDFs extracted with naive Fourier transform and withthe BG method are now qualitatively closer to the phe-nomenological distributions. The fitting reconstructionagain provides good description of the matched ITDs( χ / dof ≈ . z max /a = 8. In fact, they are also com-patible with the case of the shortest range of Ioffe time, z max /a = 4, within uncertainties. The larger range inIoffe times again decreases the uncertainty in the low- x region. For instance, at x = 0 .
2, the error in q v s is ap-proximately twice smaller with ν max ≈ . ν max ≈ .
6. An important, if somewhat obvious, effectof enhancing the ν -range is that the inverse problem inthe distribution reconstruction becomes less ill-defined.This manifests itself in the gradual convergence of all 3reconstruction methods. At z max /a = 4, the amount ofinformation in the lattice data is scarce and the rathermild additional assumption implicit in the BG method helps only little. The physically-motivated assumptionprovided by the fitting ansatz is, in turn, enough to recon-struct particularly the large- x part of the distributions.When increasing z max , there is more information from thelattice data that shapes the functional form of the dis-tributions and the mild assumption of the BG methodis enough to obtain PDFs consistent with the ones fromthe fitting ansatz. This holds in the full x -range, for allconsidered distributions. The stronger assumption con-tained in the fitting ansatz is, in turn, verified by addinglarger- ν data, proving that the estimates of large- x PDFsare robust. The access to low x < ∼ . ν max ≈ . z max /a =4 , , , ,
12. It clearly demonstrates that the naiveFourier transform is not a plausible method to recon-struct the x -dependence. All naively reconstructed distri-butions are unrobust against changing the ν -range and,at best, certain qualitative agreement is observed withphenomenological PDFs when z max is large. The BGmethod does significantly better already with intermedi-ate values of z max . The qualitative features of the PDFs1 Naive FT
Backus-Gilbert fits
FIG. 9. The dependence of thereconstructed distributions on therange of available Ioffe times,proxied by z max /a .The upper 4 plots are for q v , q v s , q , ¯ q with the naive Fourier transform,the middle 4 for the Backus-Gilbertmethod and the bottom 4 forthe fitting ansatz reconstruction. FIG. 10. The dependence of the fitting ansatz reconstructed distributions (upper left: q v , upper right: q v s , lower left: q , lowerright: ¯ q ) on the value of the strong coupling constant, α s /π . The range of Ioffe times in the reconstruction is up to ν max ≈ . z max /a = 8). are reproduced, especially at larger x . With large z max ,there is even quantitative agreement with NNPDF for awide range of x . Comparing q v and q v s , the former ismore robustly reconstructed due to the faster decay ofthe real part of matched ITDs in ν – effectively, there ismore “missing information” in the imaginary part. Thismissing data is also manifested in the irregular behaviorof the error, with PDFs at some x having artificially sup-pressed or enhanced errors. It is rather evident that themost robust way of reconstructing the PDFs is the fittingansatz approach. The numerically stronger assumptionfor regulating the inverse problem than the one in theBG method, but physically a well-motivated one, leadsto a regular behavior of PDFs with respect to the rangeof Ioffe times. Distributions obtained from all values of ν max are compatible with one another and the parameter ν max predominantly controls the uncertainties, in generaldecreasing the errors at lower values of x when data atlarger ν are included. At the same time, the large- x partsof the PDFs are little affected, since they are controlledmostly by the small- ν ITDs.As a further check of systematics, we vary the strongcoupling constant in the one-loop evolution and match-ing of ITDs. In Fig. 10, we compare our choice of this coupling, α s /π ≈ .
129 at the MS scale of 2 GeV, with α s /π = 0 .
1, equal to or close to the value used in otherstudies [64, 65]. Even though the change of α s is ratherlarge, the final distributions are not heavily affected. Thelargest effect is observed for x ≈ . − . q v s and itgets propagated also to the full distribution q . However,even this change is well within statistical uncertainties.The robustness of the PDFs with respect to α s is a conse-quence of the opposite sign of effects in the evolution andthe scheme conversion. While the evolved ITDs dependon α s in a more pronounced way, the effect of the schemeconversion brings them back towards reduced matrix el-ements and thus, the dependence of the matched ITDsand of the light-cone PDFs on α s is relatively mild. Itremains to be established whether this feature holds alsoat higher-loop orders.The final question that we want to address in thissubsection is how much the preconditioning affects theBG reconstruction. Above, we observed that the BG re-sults become increasingly consistent with the ones fromthe fitting reconstruction when ν max is increased. At ν max ≈ .
9, both methods agree within uncertaintiesfor the full x -range of all distributions, as illustrated inFig. 8. Now, we test whether this agreement can be3 FIG. 11. Comparison of the distributions (upper left: q v , upper right: q v s , lower left: q , lower right: ¯ q ) reconstructed withthe fitting ansatz approach (red) and the Backus-Gilbert method with (green) and without (blue) preconditioning. The rangeof Ioffe times in the reconstruction is up to ν max ≈ . z max /a = 8), α s /π ≈ . extended to lower ν max when preconditioning the BGmethod with the function found in the fitting ansatz ap-proach, i.e. the rescaling function p ( x ) in Eq. (17) is ofthe form (18) with parameters a, b (for q v ) or a, b, N (for q v s ) taken to be the central values found in the fits.We emphasize that this does not enforce such form ofthe distribution, but applies the BG criterion of maxi-mal stability of the solution with respect to statisticalvariance of the data to the deviation of the distributionfrom the assumed one, instead of to the full distribution.In this way, this tests the consistency of the BG assump-tion regulating the inverse problem with the assumptionmade in the fitting ansatz approach. In Fig. 11, the com-parison of the distributions from the BG method with(green band) and without (blue band) preconditioning tothe ones from the fits (red band) is given for ν max ≈ . x . However,the BG criterion for regulating the inverse problem isnot completely equivalent to the fitting ansatz assump-tion, i.e. the reconstructed function ˜ q ( x ) = q ( x ) /p ( x ) (seenotation above Eq. (17)) is not equal to 1. The statisti-cally most prominent effect of preconditioning is observed at large x – the preconditioned distributions evince nowfully smooth behavior. The agreement of PDFs fromBG precond with global fits is, however, slightly worse thanfor the fitting approach, particularly at x ≈ . q v . Aswe noted already for the standard BG method in Fig. 9,one observes irregular behavior of the errors, which areartificially suppressed for some x values and enhanced forothers. It can be interpreted that the statistical varianceof our ITDs allows only for a very restricted value of thereconstructed PDF at some x , while at other values of x significant variation is possible. Naturally, such effectsare linked to the missing data in Ioffe time. We note that ν ∈ [0 , .
2] is the minimum range of Ioffe times neededto observe consistency between BG precond and fits. Thisagain points to the fact that ν max ≈ . D. Final results with quantified systematicuncertainties
Having analyzed in detail three distribution recon-struction methods and the dependence of the results onthe range of Ioffe times and on the value of the strongcoupling constant, we are in position to present our finalPDFs. Given the robustness of the fitting ansatz recon-struction with respect to ν max and also α s , we choosethis approach as our preferred one. The central valuesof our lattice-extracted PDFs use matched ITDs in therange ν ∈ [0 , .
2] ( z max /a = 8), which fall on a universalcurve, i.e. Q ( ν, µ ) that can be obtained from differentcombinations of ( P , z ) are compatible with each otherwithin statistical uncertainties.Apart from statistical errors of thus defined PDFs, wealso add systematic uncertainties. First, we consider theuncertainty related to the range of Ioffe times to be takenin the reconstruction procedure, ∆ z max . For each distri-bution q , we use a conservative definition, that is∆ z max ( x ) = | q z max /a =12 ( x ) − q z max /a =4 ( x ) | . (20)Second, we consider the uncertainty from the choice of α s , ∆ α s ( x ):∆ α s ( x ) = | q α s /π =0 . ( x ) − q α s /π =0 . ( x ) | . (21)These two are added to the statistical error and can beconsidered as the quantified systematics of our result.In addition, there are further systematic effects in thecomputation, that cannot be quantified at the presentstage. The latter will be subject to extensive follow-upwork and will require additional numerical computationsand further theoretical developments. An extensive dis-cussion of these effects is given in Ref. [79] and below, wecomment on the most relevant points and we follow thestrategy of Ref. [49] of assuming plausible magnitudes ofthe considered effects, as percentages of ITDs values fordifferent Ioffe times. We use estimates corresponding toscenario labeled S2 in Ref. [49], considered to be the mostrealistic one.Discretization effects are an obvious source of system-atics in lattice calculations. Our results were obtained ata single value of the lattice spacing and thus, they arepotentially contaminated by these effects. To eliminatethis uncertainty, simulations at preferably at least twoadditional lattice spacings are required. Before this isdone, in a longer time perspective due to the large costof such simulations at the physical point, we assume thatcutoff effects in our present data can be up to 20%. Thisnumber is rather conservative. One argument to sup-port this claim is that larger discretization effects wouldinevitably lead to the violation of the continuum disper-sion relation, E = P + m N , where m N is the nucleonmass. Meanwhile, the dispersion relation was tested inRef. [45] and no deviations from the expected continuumbehavior were found. Additionally, related computations of moments of unpolarized PDFs by different groups (seee.g. Ref. [93]) found deviations of O (5 − O (1 − m π L > ∼
3. In our case, m π L ≈
3. However,as pointed out in Ref. [34] based on calculations in a toyscalar model, FVE in matrix element computations withspatially extended operators may be enhanced and therelevant parameter may be m N ( L − z ). Since the nucleonmass is much larger than the pion mass, it would effec-tively not lead to any enhancement of FVE. The worstplausible case is if in QCD the parameter controlling FVEbecomes m π ( L − z ). However, even then, the values of z that are actually used should not lead to severe FVE.This was confirmed for the z -dependent renormalizationfunctions used in the non-perturbative renormalizationof quasi-PDFs in Ref. [45]. Overall, to remain conserva-tive, we allow for 5% FVE in our hypothetical systematicerror budget.Further systematic effects can result from contamina-tion of the signal for the nucleon by excited states withthe same quantum numbers. For this uncertainty, inves-tigated in great detail for the matrix elements used inthis work, we rely on the conclusion of Ref. [45], whereexcited states suppression was found within statisticalerrors of the results. Thus, we take this kind of system-atics to be 10% for all ITDs, i.e. slightly larger than theattained statistical precision.The perturbative ingredient of our computation is sub-ject to truncation effects. These need to be investigatedby calculating at least the two-loop matching. Until thisis carried out, the magnitude of this systematic effectis unknown. As we have demonstrated above, chang-ing the α s value for the one-loop formula does not leadto large changes of the PDFs. However, the neglectedhigher-order effects may still be sizable. The matching tolight-cone ITDs, i.e. the factorization of the pseudo-ITDinto its light-cone counterpart and a perturbative coeffi-cient is also subject to higher-twist effects of O ( z Λ ).These effects, again, need to be further investigated withdedicated calculations, but they are not expected to beoverwhelming with the values of z that are included inthe analysis and in view of the rather mild dependenceon z max that we find in this work. Overall, for pointsdiscussed in this paragraph, we conservatively attribute20% as their potential size.Our final PDFs are shown in Fig. 12 and comparedto global fits of NNPDF [92]. We show 3 kinds of er-ror bands. The purple one (the most narrow) is exclu-sively the statistical error of our results. The systematicuncertainties, discussed above, enter in the blue band(quantified systematics from varying the range of Ioffetimes and the value of the strong coupling) and in the5 FIG. 12. Final unpolarized PDFs extracted from the lattice (fitting ansatz reconstruction) and compared to global fits ofNNPDF [92] (solid black line and dark grey band). Shown distributions: valence ( q v ; upper left), valence + 2 sea ( q v s = q v +2¯ q ;upper right), full ( q = q v + ¯ q ; lower left) and sea (¯ q ; lower right). The range of Ioffe times in the reconstruction is up to ν max ≈ . z max /a = 8), α s /π ≈ . ν max and α s (blue) and the totalerror additionally with estimated uncertainties related to cutoff effects, FVE, excited states contamination, truncation andhigher-twist effects (cyan) – see text for more details. cyan band (conservatively estimated errors from cutoff ef-fects, FVE, excited states contamination, truncation andhigher-twist effects). The total uncertainty combines allthe separate sources thereof in quadrature.For all distributions, we find very good agreementwith the corresponding phenomenological curve alreadywithin statistical errors, while the total error accountsfor the remaining small discrepancies in certain regionsof x in q v and q (around x = 0 . x < ∼ .
05 for q v ). This gives confidence in the estimates of unquan-tified systematics, but we emphasize that much work isneeded to properly quantify these effects.The agreement of our final PDFs with NNPDF is strik-ing and shows that a lattice extraction of the full x -dependence of PDFs is feasible. It allows us also todraw conclusions about the reliability of such an extrac-tion in different regions of x . It is clear that the large- x part ( x > ∼ .
6) is reconstructed very robustly for alldistributions. Since this region is dominated by small- ν ITDs, it is basically insensitive to the range of Ioffe timesand moreover, one can argue that some other sources of systematics are expected to be small. For instance, O ( z Λ ) higher-twist effects are probably negligible,there are definitely no enhanced FVE of the type dis-cussed in Ref. [34] and cutoff effects are likely small,since the calculated small- z matrix elements do not differmuch from the local ( z = 0) ones, for which automatic O ( a )-improvement holds in our setup. In the intermedi-ate range of x , x ≈ . − .
6, the uncertainties tend toincrease and the one from varying the Ioffe time rangebecomes non-negligible, even if it is still subleading. Thereconstructed PDFs are less constrained due to part ofthe Ioffe time dependence of ITDs missing. The latter be-comes especially important in the low- x region ( x < ∼ . V. SUMMARY
In this work, we used the pseudo-distribution approachto calculate the x -dependence of the unpolarized par-ton distribution functions of the nucleon. The methodrelies on computations of spatial correlations betweenboosted nucleon states. The resulting matrix elementsare then used to form appropriate ratios that cancel theexisting logarithmic and power-like divergences, definingpseudo-distributions in Ioffe time (Ioffe time distribu-tions or ITDs). Pseudo-ITDs are then matched to theirlight-cone counterparts, objects containing physical in-formation and Fourier-conjugate to PDFs. The step oftranslating the light-cone ITDs to PDFs is highly non-trivial due to the difficulty of obtaining the full Ioffe-timedependence of the former. Thus, it is subject to an in-verse problem and advanced reconstruction methods areneeded instead of a simple Fourier transform.We used a robust lattice setup of maximally twistedmass fermions to compute the bare matrix elements atthe physical pion mass and with a relatively fine lat-tice spacing and large lattice volume that should leadto at most modest discretization and finite volume ef-fects. Having computed pseudo-ITDs, we performed thematching procedure to go to the light cone and we re-constructed the x -distributions using 3 methods. As ex-pected, the naive Fourier transform does not lead to ro-bust results. The Backus-Gilbert methods offers one so-lution to the inverse problem by assuming that the re-constructed distribution should have minimal variancewith respect to the statistical variation of the data. Wefound that with large range of the Ioffe-time depen-dence missing, this criterion is not good enough to re-construct PDFs. However, with this range being ex-tended, the method gives results convergent with theones from the third reconstruction approach. The lat-ter assumes a functional form of the light-cone PDFs, asdone in phenomenological analyses to extract PDFs fromglobal fits. This assumption, similarly as the one in theBackus-Gilbert method, regulates the inverse problemand leads, in practice, to well-behaved and robust PDFs.We checked this robustness by investigating the depen-dence of the reconstructed distributions with respect tothe range of included Ioffe times and also the value ofthe strong coupling constant used in the matching. Wefound that the large- x region of PDFs is insensitive tothis range, while at smaller x , ITDs at large Ioffe timelead to variations of the final PDFs only within statis-tical uncertainties. Obviously, the missing data at theselarge Ioffe times increase the error of PDF estimates atlow- x . To obtain more precise results there, ITDs needto be computed at larger nucleon boosts. This presents apractical problem for the lattice, since the signal-to-noiseratio in the computation of matrix elements is quickly de-caying with each additional unit of momentum. Thus, itis natural that it will be difficult for the lattice to reli-ably extract the very low- x region. Nevertheless, a largepart of the x -dependence is possible to determine on the lattice with relatively modest computational resources.Despite the optimistic results obtained here, it needsto be bore in mind that fully robust lattice-extractedPDFs need to have all relevant sources of systematic ef-fects quantified. This will need a significant load of workin the next few years, involving additional computationson the lattice as well as theoretical developments. Somesources of systematics can be quantified in a straight-forward manner by repeating the procedure used in thiswork for additional ensembles of gauge field configura-tions, in particular at finer lattice spacings and largervolumes. This will allow us to address, respectively, dis-cretization and finite volume effects. For other kinds ofeffects, proper strategies of addressing them need to bedevised. Truncation effects can only be properly quan-tified after derivation of higher-loop matching. Higher-twist effects, in turn, can, in principle, be accessed nu-merically by computing proper matrix elements of morecomplicated operators, but analytical insight can be in-valuable as well.Thus, the field of extracting x -dependent PDFs, as wellas other partonic functions, still needs a lot of progress.However, it is clearly reassuring that the quality of ourpresent results is already very satisfactory. Only a fewyears ago, lattice data for PDFs were limited to only thelowest two or three moments, without realistic perspec-tives of reconstructing the x -dependence. The progressinduced by the seminal paper of Ji proposing how toextract the latter and subsequent alternative proposals,in particular the one we used in this work, change theprospects of this field to a huge extent. As demonstratedin this work, the PDFs can indeed be extracted directlyfrom first principles , i.e. from the QCD Lagrangian andalready now this can be done with both qualitative andquantitative agreement with global fits. The latter re-quires additional work to be fully established and at thisstage, we resorted to plausible hypotheses about the sizeof some systematic effects.Given the success of this program for unpolarizedPDFs, obvious directions for the future, apart from thediscussed analysis of systematics, is to extend the workto polarized PDFs and other kinds of structure functions,in particular generalized parton distributions (GPDs)and transverse-momentum-dependent parton distribu-tions (TMDs), as well as to singlet distributions. Allof these directions are challenging, but offer an unprece-dented opportunity of having crucial insights for the par-tonic structure of the nucleon, relevant for its deeperunderstanding both at the theoretical and experimentallevel. ACKNOWLEDGMENTS
K.C. thanks Savvas Zafeiropoulos for many interestingdiscussions about the pseudo-distribution approach atdifferent stages of evolution of this method. M.B., K.C.and A.S. are supported by the National Science Centre7(Poland) grant SONATA BIS no. 2016/22/E/ST2/00013.M.C. acknowledges financial support by the U.S. Depart-ment of Energy, Office of Nuclear Physics, Early CareerAward under Grant No. de-sc0020405. This research wassupported in part by PLGrid Infrastructure (Prometheussupercomputer at AGH Cyfronet in Cracow). Com- putations were also partially performed at the PoznanSupercomputing and Networking Center (Eagle super-computer), the Interdisciplinary Centre for Mathemat-ical and Computational Modelling of the Warsaw Uni-versity (Okeanos supercomputer) and at the AcademicComputer Centre in Gda´nsk (Tryton supercomputer). [1] N. A. of Sciences Engineering Medicine,
An Assessmentof U.S.-Based Electron-Ion Collider Science (2018).[2] A. Abdel-Rehim et al. , Nucleon and pion structure withlattice QCD simulations at physical value of the pionmass, Phys. Rev.
D92 , 114513 (2015), [Erratum: Phys.Rev.D93,no.3,039904(2016)], arXiv:1507.04936 [hep-lat].[3] M. Oehm, C. Alexandrou, M. Constantinou, K. Jansen,G. Koutsou, B. Kostrzewa, F. Steffens, C. Urbach, andS. Zafeiropoulos, (cid:104) x (cid:105) and (cid:104) x (cid:105) of the pion PDF from lat-tice QCD with N f = 2 + 1 + 1 dynamical quark flavors,Phys. Rev. D99 , 014508 (2019), arXiv:1810.09743 [hep-lat].[4] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finken-rath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, andA. Vaquero Aviles-Casco, Proton and neutron electro-magnetic form factors from lattice QCD, Phys. Rev.
D100 , 014509 (2019), arXiv:1812.10311 [hep-lat].[5] C. Alexandrou et al. , Moments of nucleon generalizedparton distributions from lattice QCD simulations atphysical pion mass, Phys. Rev.
D101 , 034519 (2020),arXiv:1908.10706 [hep-lat].[6] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finken-rath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, andA. Vaquero Aviles-Casco, The nucleon axial, tensor andscalar charges and σ -terms in lattice QCD, (2019),arXiv:1909.00485 [hep-lat].[7] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finken-rath, K. Hadjiyiannakou, K. Jansen, and G. Koutsou,Nucleon strange electromagnetic form factors, Phys. Rev. D101 , 031501 (2020), arXiv:1909.10744 [hep-lat].[8] X. Ji, Parton Physics on a Euclidean Lattice,Phys.Rev.Lett. , 262002 (2013), arXiv:1306.1539[hep-ph].[9] X. Ji, Parton Physics from Large-Momentum EffectiveField Theory, Sci. China Phys. Mech. Astron. , 1407(2014), arXiv:1404.6680 [hep-ph].[10] X. Xiong, X. Ji, J.-H. Zhang, and Y. Zhao, One-loop matching for parton distributions: Nonsinglet case,Phys.Rev. D90 , 014051 (2014), arXiv:1310.7471 [hep-ph].[11] H.-W. Lin, J.-W. Chen, S. D. Cohen, and X. Ji, FlavorStructure of the Nucleon Sea from Lattice QCD, Phys.Rev.
D91 , 054510 (2015), arXiv:1402.1462 [hep-ph].[12] L. Gamberg, Z.-B. Kang, I. Vitev, and H. Xing,Quasi-parton distribution functions: a study in the di-quark spectator model, Phys. Lett.
B743 , 112 (2015),arXiv:1412.3401 [hep-ph].[13] C. Alexandrou, K. Cichy, V. Drach, E. Garcia-Ramos,K. Hadjiyiannakou, K. Jansen, F. Steffens, and C. Wiese,Lattice calculation of parton distributions, Phys. Rev.
D92 , 014502 (2015), arXiv:1504.07455 [hep-lat].[14] I. Vitev, L. Gamberg, Z. Kang, and H. Xing, A Studyof Quasi-parton Distribution Functions in the Diquark Spectator Model,
Proceedings, QCD Evolution Workshop(QCD 2015): Newport News, VA, USA, May 26-30,2015 , PoS
QCDEV2015 , 045 (2015), arXiv:1511.05242[hep-ph].[15] Y. Jia and X. Xiong, Quasidistribution amplitude ofheavy quarkonia, Phys. Rev.
D94 , 094005 (2016),arXiv:1511.04430 [hep-ph].[16] J.-W. Chen, S. D. Cohen, X. Ji, H.-W. Lin, and J.-H.Zhang, Nucleon Helicity and Transversity Parton Dis-tributions from Lattice QCD, Nucl. Phys.
B911 , 246(2016), arXiv:1603.06664 [hep-ph].[17] J.-W. Chen, X. Ji, and J.-H. Zhang, Improved quasiparton distribution through Wilson line renormalization,Nucl. Phys.
B915 , 1 (2017), arXiv:1609.08102 [hep-ph].[18] C. Alexandrou, K. Cichy, M. Constantinou, K. Had-jiyiannakou, K. Jansen, F. Steffens, and C. Wiese, Up-dated Lattice Results for Parton Distributions, Phys.Rev.
D96 , 014513 (2017), arXiv:1610.03689 [hep-lat].[19] A. Bacchetta, M. Radici, B. Pasquini, and X. Xiong,Reconstructing parton densities at large fractional mo-menta, Phys. Rev.
D95 , 014036 (2017), arXiv:1608.07638[hep-ph].[20] R. A. Brice˜no, M. T. Hansen, and C. J. Monahan, Roleof the Euclidean signature in lattice calculations of qua-sidistributions and other nonlocal matrix elements, Phys.Rev.
D96 , 014502 (2017), arXiv:1703.06072 [hep-lat].[21] T. Ishikawa, Y.-Q. Ma, J.-W. Qiu, and S. Yoshida,Renormalizability of quasiparton distribution functions,Phys. Rev.
D96 , 094019 (2017), arXiv:1707.03107 [hep-ph].[22] X. Ji, J.-H. Zhang, and Y. Zhao, Renormalization inLarge Momentum Effective Theory of Parton Physics,Phys. Rev. Lett. , 112001 (2018), arXiv:1706.08962[hep-ph].[23] M. Constantinou and H. Panagopoulos, Perturbativerenormalization of quasi-parton distribution functions,Phys. Rev.
D96 , 054506 (2017), arXiv:1705.11193 [hep-lat].[24] C. Alexandrou, K. Cichy, M. Constantinou, K. Had-jiyiannakou, K. Jansen, H. Panagopoulos, and F. Stef-fens, A complete non-perturbative renormalization pre-scription for quasi-PDFs, Nucl. Phys.
B923 , 394 (2017),arXiv:1706.00265 [hep-lat].[25] X. Ji, J.-H. Zhang, and Y. Zhao, More On Large-Momentum Effective Theory Approach to PartonPhysics, Nucl. Phys.
B924 , 366 (2017), arXiv:1706.07416[hep-ph].[26] W. Wang, S. Zhao, and R. Zhu, Gluon quasidistribu-tion function at one loop, Eur. Phys. J.
C78 , 147 (2018),arXiv:1708.02458 [hep-ph].[27] J. Green, K. Jansen, and F. Steffens, NonperturbativeRenormalization of Nonlocal Quark Bilinears for Par-ton Quasidistribution Functions on the Lattice Using an Auxiliary Field, Phys. Rev. Lett. , 022004 (2018),arXiv:1707.07152 [hep-lat].[28] I. W. Stewart and Y. Zhao, Matching the quasipartondistribution in a momentum subtraction scheme, Phys.Rev. D , 054512 (2018), arXiv:1709.04933 [hep-ph].[29] W. Broniowski and E. Ruiz Arriola, Nonperturbativepartonic quasidistributions of the pion from chiral quarkmodels, Phys. Lett. B773 , 385 (2017), arXiv:1707.09588[hep-ph].[30] W. Broniowski and E. Ruiz Arriola, Partonic qua-sidistributions of the proton and pion from transverse-momentum distributions, Phys. Rev.
D97 , 034031(2018), arXiv:1711.03377 [hep-ph].[31] T. Izubuchi, X. Ji, L. Jin, I. W. Stewart, and Y. Zhao,Factorization Theorem Relating Euclidean and Light-Cone Parton Distributions, Phys. Rev.
D98 , 056004(2018), arXiv:1801.03917 [hep-ph].[32] C. Alexandrou, K. Cichy, M. Constantinou, K. Jansen,A. Scapellato, and F. Steffens, Light-Cone Parton Dis-tribution Functions from Lattice QCD, Phys. Rev. Lett. , 112001 (2018), arXiv:1803.02685 [hep-lat].[33] J.-H. Zhang, J.-W. Chen, L. Jin, H.-W. Lin, A. Sch¨afer,and Y. Zhao, First direct lattice-QCD calculation of the x -dependence of the pion parton distribution function,Phys. Rev. D100 , 034505 (2019), arXiv:1804.01483 [hep-lat].[34] R. A. Brice˜no, J. V. Guerrero, M. T. Hansen, andC. J. Monahan, Finite-volume effects due to spatiallynonlocal operators, Phys. Rev.
D98 , 014511 (2018),arXiv:1805.01034 [hep-lat].[35] G. Spanoudes and H. Panagopoulos, Renormaliza-tion of Wilson-line operators in the presence ofnonzero quark masses, Phys. Rev.
D98 , 014509 (2018),arXiv:1805.01164 [hep-lat].[36] A. Radyushkin, Structure of parton quasi-distributionsand their moments, Phys. Lett. B , 380 (2019),arXiv:1807.07509 [hep-ph].[37] J. Karpie, K. Orginos, and S. Zafeiropoulos, Moments ofIoffe time parton distribution functions from non-localmatrix elements, JHEP , 178, arXiv:1807.10933 [hep-lat].[38] Y.-S. Liu et al. (Lattice Parton), Unpolarized isovectorquark distribution function from lattice QCD: A system-atic analysis of renormalization and matching, Phys. Rev.D , 034020 (2020), arXiv:1807.06566 [hep-lat].[39] H.-W. Lin, J.-W. Chen, X. Ji, L. Jin, R. Li, Y.-S. Liu, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Proton Isovector He-licity Distribution on the Lattice at Physical Pion Mass,Phys. Rev. Lett. , 242003 (2018), arXiv:1807.07431[hep-lat].[40] Y. Jia, S. Liang, X. Xiong, and R. Yu, Partonic quasidis-tributions in two-dimensional QCD, Phys. Rev. D98 ,054011 (2018), arXiv:1804.04644 [hep-th].[41] S. Bhattacharya, C. Cocuzza, and A. Metz, Generalizedquasi parton distributions in a diquark spectator model,Phys. Lett. B , 453 (2019), arXiv:1808.01437 [hep-ph].[42] C. Alexandrou, K. Cichy, M. Constantinou, K. Jansen,A. Scapellato, and F. Steffens, Transversity parton dis-tribution functions from lattice QCD, Phys. Rev. D ,091503 (2018), arXiv:1807.00232 [hep-lat].[43] V. M. Braun, A. Vladimirov, and J.-H. Zhang, Powercorrections and renormalons in parton quasidistributions,Phys. Rev. D99 , 014013 (2019), arXiv:1810.00048 [hep- ph].[44] J. Karpie, K. Orginos, A. Rothkopf, and S. Zafeiropoulos,Reconstructing parton distribution functions from Ioffetime data: from Bayesian methods to Neural Networks,JHEP , 057, arXiv:1901.05408 [hep-lat].[45] C. Alexandrou, K. Cichy, M. Constantinou, K. Had-jiyiannakou, K. Jansen, A. Scapellato, and F. Steffens,Systematic uncertainties in parton distribution functionsfrom lattice QCD simulations at the physical point, Phys.Rev. D99 , 114504 (2019), arXiv:1902.00587 [hep-lat].[46] S. Bhattacharya, C. Cocuzza, and A. Metz, Exploringtwist-2 GPDs through quasi-distributions in a diquarkspectator model, (2019), arXiv:1903.05721 [hep-ph].[47] J.-W. Chen, H.-W. Lin, and J.-H. Zhang, Piongeneralized parton distribution from latticeQCD 10.1016/j.nuclphysb.2020.114940 (2019),arXiv:1904.12376 [hep-lat].[48] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukher-jee, P. Petreczky, C. Shugert, and S. Syritsyn, Valenceparton distribution function of pion from fine lattice,Phys. Rev.
D100 , 034516 (2019), arXiv:1905.06349 [hep-lat].[49] K. Cichy, L. Del Debbio, and T. Giani, Parton distribu-tions from lattice data: the nonsinglet case, JHEP ,137, arXiv:1907.06037 [hep-ph].[50] H.-D. Son, A. Tandogan, and M. V. Polyakov, Nucleonquasi-Parton Distributions in the large N c limit, (2019),arXiv:1911.01955 [hep-ph].[51] Y. Chai et al. , Parton distribution functions of ∆ + onthe lattice, (2020), arXiv:2002.12044 [hep-lat].[52] J. R. Green, K. Jansen, and F. Steffens, Improvement,generalization, and scheme conversion of Wilson-line op-erators on the lattice in the auxiliary field approach,Phys. Rev. D , 074509 (2020), arXiv:2002.09408 [hep-lat].[53] X. Ji, Fundamental Properties of the Proton in Light-Front Zero Modes, (2020), arXiv:2003.04478 [hep-ph].[54] S. Bhattacharya, K. Cichy, M. Constantinou, A. Metz,A. Scapellato, and F. Steffens, New insights on protonstructure from lattice QCD: the twist-3 parton distribu-tion function g T ( x ), (2020), arXiv:2004.04130 [hep-lat].[55] A. Radyushkin, Nonperturbative Evolution of PartonQuasi-Distributions, Phys. Lett. B767 , 314 (2017),arXiv:1612.05170 [hep-ph].[56] A. V. Radyushkin, Quasi-parton distribution func-tions, momentum distributions, and pseudo-parton dis-tribution functions, Phys. Rev.
D96 , 034025 (2017),arXiv:1705.01488 [hep-ph].[57] A. V. Radyushkin, Quark pseudodistributions atshort distances, Phys. Lett.
B781 , 433 (2018),arXiv:1710.08813 [hep-ph].[58] A. Radyushkin, Quasi-PDFs and pseudo-PDFs,
Proceed-ings, QCD Evolution Workshop (QCD 2017): NewportNews, VA, USA, May 22-26, 2017 , PoS
QCDEV2017 ,021 (2017), arXiv:1711.06031 [hep-ph].[59] A. Radyushkin, One-loop evolution of parton pseudo-distribution functions on the lattice, Phys. Rev.
D98 ,014019 (2018), arXiv:1801.02427 [hep-ph].[60] A. V. Radyushkin, Generalized parton distributions andpseudodistributions, Phys. Rev. D , 116011 (2019),arXiv:1909.08474 [hep-ph].[61] A. V. Radyushkin, Theory and applications of par-ton pseudodistributions, International Journal of ModernPhysics A , 2030002 (2020), arXiv:1912.04244 [hep- ph].[62] B. L. Ioffe, Space-time picture of photon and neutrinoscattering and electroproduction cross-section asymp-totics, Phys. Lett. , 123 (1969).[63] J.-H. Zhang, J.-W. Chen, and C. Monahan, Parton dis-tribution functions from reduced Ioffe-time distributions,Phys. Rev. D97 , 074508 (2018), arXiv:1801.03023 [hep-ph].[64] K. Orginos, A. Radyushkin, J. Karpie, and S. Zafeiropou-los, Lattice QCD exploration of parton pseudo-distribution functions, Phys. Rev.
D96 , 094503 (2017),arXiv:1706.05373 [hep-ph].[65] B. Jo´o, J. Karpie, K. Orginos, A. Radyushkin,D. Richards, and S. Zafeiropoulos, Parton DistributionFunctions from Ioffe time pseudo-distributions, JHEP , 081, arXiv:1908.09771 [hep-lat].[66] B. Jo´o, J. Karpie, K. Orginos, A. V. Radyushkin,D. G. Richards, R. S. Sufian, and S. Zafeiropoulos,Pion valence structure from Ioffe-time parton pseudodis-tribution functions, Phys. Rev. D , 114512 (2019),arXiv:1909.08517 [hep-lat].[67] B. Jo´o, J. Karpie, K. Orginos, A. V. Radyushkin,D. G. Richards, and S. Zafeiropoulos, Parton Distri-bution Functions from Ioffe time pseudo-distributionsfrom lattice caclulations; approaching the physical point,(2020), arXiv:2004.01687 [hep-lat].[68] U. Aglietti, M. Ciuchini, G. Corbo, E. Franco, G. Mar-tinelli, and L. Silvestrini, Model independent determi-nation of the light cone wave functions for exclusiveprocesses, Phys. Lett. B441 , 371 (1998), arXiv:hep-ph/9806277 [hep-ph].[69] W. Detmold and C. J. D. Lin, Deep-inelastic scatter-ing and the operator product expansion in lattice QCD,Phys. Rev.
D73 , 014501 (2006), arXiv:hep-lat/0507007[hep-lat].[70] V. Braun and D. Mueller, Exclusive processes in positionspace and the pion distribution amplitude, Eur. Phys. J.
C55 , 349 (2008), arXiv:0709.1348 [hep-ph].[71] Y.-Q. Ma and J.-W. Qiu, Extracting Parton DistributionFunctions from Lattice QCD Calculations, Phys. Rev.
D98 , 074021 (2018), arXiv:1404.6860 [hep-ph].[72] Y.-Q. Ma and J.-W. Qiu, Exploring Partonic Structureof Hadrons Using ab initio Lattice QCD Calculations,Phys. Rev. Lett. , 022003 (2018), arXiv:1709.03018[hep-ph].[73] K.-F. Liu and S.-J. Dong, Origin of difference betweenanti-d and anti-u partons in the nucleon, Phys. Rev. Lett. , 1790 (1994), arXiv:hep-ph/9306299 [hep-ph].[74] A. J. Chambers, R. Horsley, Y. Nakamura, H. Perlt,P. E. L. Rakow, G. Schierholz, A. Schiller, K. Somfleth,R. D. Young, and J. M. Zanotti, Nucleon Structure Func-tions from Operator Product Expansion on the Lattice,Phys. Rev. Lett. , 242001 (2017), arXiv:1703.01153[hep-lat].[75] G. S. Bali, V. M. Braun, B. Gl¨aßle, M. G¨ockeler,M. Gruber, F. Hutzler, P. Korcyl, A. Sch¨afer, P. Wein,and J.-H. Zhang, Pion distribution amplitude from Eu-clidean correlation functions: Exploring universality andhigher-twist effects, Phys. Rev. D98 , 094507 (2018),arXiv:1807.06671 [hep-lat].[76] W. Detmold, I. Kanamori, C. D. Lin, S. Mondal, andY. Zhao, Moments of pion distribution amplitude usingoperator product expansion on the lattice, PoS
LAT-TICE2018 , 106 (2018), arXiv:1810.12194 [hep-lat]. [77] R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu,and D. G. Richards, Pion Valence Quark Distributionfrom Matrix Element Calculated in Lattice QCD, Phys.Rev.
D99 , 074507 (2019), arXiv:1901.03921 [hep-lat].[78] J. Liang, T. Draper, K.-F. Liu, A. Rothkopf, and Y.-B. Yang (XQCD), Towards the nucleon hadronic tensorfrom lattice QCD, (2019), arXiv:1906.05312 [hep-ph].[79] K. Cichy and M. Constantinou, A guide to light-conePDFs from Lattice QCD: an overview of approaches,techniques and results, Adv. High Energy Phys. ,3036904 (2019), arXiv:1811.07248 [hep-lat].[80] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Large-Momentum Effective Theory, (2020), arXiv:2004.03543[hep-ph].[81] G. Backus and F. Gilbert, The resolving power of grossearth data, Geophysical Journal International , 169(1968).[82] A. N. Tikhonov, Solution of incorrectly formulated prob-lems and the regularization method, Soviet Math. Dokl. , 1035 (1963).[83] M. V. Ulybyshev, C. Winterowd, and S. Zafeiropoulos,Direct detection of metal-insulator phase transitions us-ing the modified Backus-Gilbert method, EPJ Web Conf. , 03008 (2018), arXiv:1710.06675 [hep-lat].[84] M. Ulybyshev, C. Winterowd, and S. Zafeiropoulos, Col-lective charge excitations and the metal-insulator tran-sition in the square lattice Hubbard-Coulomb model,Phys. Rev. B , 205115 (2017), arXiv:1707.04212 [cond-mat.str-el].[85] R. Frezzotti, P. A. Grassi, S. Sint, and P. Weisz (Alpha),Lattice QCD with a chirally twisted mass term, JHEP , 058, arXiv:hep-lat/0101001 [hep-lat].[86] R. Frezzotti and G. C. Rossi, Chirally improving Wil-son fermions. 1. O(a) improvement, JHEP , 007,arXiv:hep-lat/0306014 [hep-lat].[87] B. Sheikholeslami and R. Wohlert, Improved ContinuumLimit Lattice Action for QCD with Wilson Fermions,Nucl. Phys. B259 , 572 (1985).[88] A. Abdel-Rehim et al. (ETM), First physics results atthe physical pion mass from N f = 2 Wilson twistedmass fermions at maximal twist, Phys. Rev. D95 , 094515(2017), arXiv:1507.05068 [hep-lat].[89] Y. Iwasaki, Renormalization Group Analysis of Lat-tice Theories and Improved Lattice Action. II. Four-dimensional non-Abelian SU(N) gauge model, (1983),arXiv:1111.7054 [hep-lat].[90] C. Alexandrou and C. Kallidonis, Low-lying baryonmasses using N f = 2 twisted mass clover-improvedfermions directly at the physical pion mass, Phys. Rev. D96 , 034511 (2017), arXiv:1704.02647 [hep-lat].[91] C. Alexandrou, M. Constantinou, and H. Panagopoulos(ETM), Renormalization functions for Nf=2 and Nf=4twisted mass fermions, Phys. Rev.
D95 , 034505 (2017),arXiv:1509.00213 [hep-lat].[92] R. D. Ball et al. (NNPDF), Parton distributions fromhigh-precision collider data, Eur. Phys. J.
C77 , 663(2017), arXiv:1706.00428 [hep-ph].[93] M. Constantinou, Hadron Structure,
Proceedings, 32ndInternational Symposium on Lattice Field Theory (Lat-tice 2014): Brookhaven, NY, USA, June 23-28, 2014 ,PoS