Bayesian-Wilson coefficients in lattice QCD computations of valence PDFs and GPDs

We propose an analysis method for the leading-twist operator product expansion based lattice QCD determinations of the valence parton distribution function (PDF). In the first step, we determine the confidence-intervals of the leading-twist $\overline{\mathrm{MS}}$ Wilson coefficients, $C_n(\mu^2 z^2)$, of the equal-time bilocal quark bilinear, given the lattice QCD matrix element of Ioffe-time distribution for a particular hadron $H$ as well as the prior knowledge of the valence PDF, $f(x,\mu)$ of the hadron $H$ determined via global fit from the experimental data. In the next step, we apply the numerically estimated $C_n$ in the lattice QCD determinations of the valence PDFs of other hadrons, and for the zero-skewness generalized parton distribution (GPD) of the same hadron $H$ at non-zero momentum transfers. Our proposal still assumes the dominance of leading-twist terms, but it offers a pragmatic alternative to the usage of perturbative Wilson coefficients and their associated higher-loop uncertainties such as the effect of all-order logarithms at larger sub-Fermi quark-antiquark separations $z$.


I. INTRODUCTION
The progress in determining the x-dependent hadron structures, such as the parton distribution functions (PDFs) and the generalized parton distribution functions (GPDs) [1][2][3] (for a review, see [4]), has been quite rapid in the recent years, particularly owing to the perturbative matching frameworks such as the quasi-PDF [5,6], pseudo-PDF [7,8], current-current correlators [9] and the good lattice cross-sections approach [10,11]. There are other ideas that are being used (for example, [12][13][14][15][16]), but we would restrict ourselves to the set of perturbative matching approaches for calculating PDFs in this paper. Roughly, the recurring idea in all the variants of the perturbative matching approaches is the factorization of certain equal-time boosted hadron matrix elements with a spatial operator separation z (or the Fourier transform thereof into x-space) into a matching kernel and the structure functions, such as the PDF as a specific case. The usefulness of this factorization is that the matching kernel can be perturbatively evaluated and hence universal for all hadrons. We refer the reader to the recent reviews [17][18][19][20] on this broad topic for further details on the methodology and for the recent works in this direction.
Without any loss of generality, in this paper, we will specifically consider the case of equal-time multiplicatively renormalizable [21,22] correlator constructed from the quasi-PDF [5] matrix element 2E(P z )M(ω, z 2 ) = P z |ψ z γ t W z,0 ψ 0 |P z , of quark and antiquark separated by a spatial distance z, also known as pseudo-Ioffe-time distribution [7] (abbreviated as pseudo-ITD or just as ITD), and evaluated in an on-shell hadron with P = (E, P z ) moving with a spatial momentum P z . The variable ω = zP z is also known as the Ioffe-time [23,24]. However, the rest of the paper can be equally carried over to any of the operators (e.g., [25,26]) under the umbrella of the good lattice cross-sections as well. The pertinent detail for this paper is that the the underlying factorization of M(ω, z 2 ), can be equivalently rewritten [9,27] as a leading twist operator product expansion (OPE), (2) valid up to higher-twist corrections, usually denoted as O(Λ 2 QCD z 2 ) terms. Occupying the central role in the above expression are the Wilson coefficients C n (z 2 µ 2 ) that match the renormalized Euclidean matrix element in some scheme to the MS PDF via the moments, and can be perturbatively computed at the hard scale set by 1/z. In a typical lattice computation, the analytic form of C n , C n (µ 2 z 2 ) =c 0 (n)+c 1 (n) ln µ 2 z 2 +c 2 (n) ln 2 (µ 2 z 2 )+. . . , (3) in terms of the coefficientsc of the powers of the singular logarithms ln µ 2 z 2 , up to certain order in α s is provided as an input, and in return, the PDF f (x, µ) and the Mellin moments, x n = 1 0 x n f (x)dx, are computed as outputs 1 .
In the rest of the paper, we will consider the nonsinglet, valence parton distribution functions. Unlike the 1 For the flavor non-singlet PDF case that we will specifically consider in this paper, the notation for the Mellin moments as given in the community white paper [28] are x n u + −d + for the odd values of n, and x n u − −d − for the even values of n that corresponds to the valence PDF moments. We will simply refer to these non-singlet moments, specifically the valence moments, as x n throughout this paper. singlet distributions, for the valence case, xf (x) vanishes in the small-x limit, and hence, enable us to neglect the associated systematic errors in the global fit determinations and in their first few Mellin moments. This lets us focus only on sources of uncertainty in the lattice computations of PDF. In the valence sector, the perturbative matching approach has been successfully applied to the lattice computations of PDF of the pion [25,26,[29][30][31] and the proton [8,[32][33][34][35][36][37][38]. In these calculations, the resulting reconstructed PDFs capture more or less the features of the phenomenologically determined PDFs (which we henceforth simply refer to as "pheno-PDFs"). By avoiding the issues of inversion problem to obtain xdependent PDF, a direct comparison of the lattice ITD data itself with the corresponding expectation from the global analysis fits, obtained via perturbative matching, have only been qualitatively close enough (for example, Ref. [34] being one such comparison with NNPDF result [39] for the proton at physical point, and Ref. [30] for comparison of JAM18 result [40] with 300 MeV pion). In particular, the deviation between the lattice ITD and the pheno-ITD increases with ω or z as discussed in [41]. These presumably point to the presence of yet uncontrolled systematic effect(s). Some of the reasons for this could be: 1. The effect of heavier than physical pion masses used in most of the calculations, but as we mentioned above, a recent computation [34] of the isovector nucleon PDF at the physical point still shows visible discrepancy with the ITD constructed from NNPDF global analysis [39]. The weak pion mass dependences seen in [26,32] in the range, m π ∈ [172, 415] MeV, suggests this might not be the dominant contribution to the observed discrepancies even for computations at moderately heavier pion masses.
2. There could be lattice spacing corrections to pseudo-ITDs. For the case of pseudo-ITD renormalized in the ratio scheme, the lattice spacing effects were seen [30] to mostly affect regions of M(ω, z 2 ) for z = O(a). Such lattice spacing effect at small z was also seen for the nucleon recently [42]. Therefore, at least for renormalization methods like the ratio, the short-distance lattice spacing effects are less important and can be ignored provided one skips first few values of z. Hence, we do not consider this in this paper.
3. A more troublesome reason could be the presence of yet undetected large higher-twist contributions in the range of separations, z ∈ [0.04, 1] fm that are accessible in the computations today, but incorrectly taken as part of leading-twist contribution. However, a recent investigation [30] of the higher-twist effect that dominate in the P z = 0 RI-MOM matrix element for the pion suggests that the higher-twist effect could be a minor, ≈ 40 MeV effect that is atypically small compared to the QCD scale, Λ QCD ≈ 300 MeV. By using ratio method (with ratio taken with respect to zero momentum as well as with non-zero momentum bare matrix element) for renormalization instead of the RI-MOM scheme [43], it is reasonable to expect these leading z 2 higher-twist contributions to get approximately canceled leaving behind only a residual higher-twist contamination of the form z 2 ω 2 [44].
In this work, we take an optimistic stand on this matter and ask what if the remaining uncertainties in the leading-twist approach are predominantly unaccounted higher-loop corrections, and how to take care of our limited knowledge of the perturbative convergence of the Wilson coefficients in a probabilistic manner? This is the motivation for this paper.
There have been previous works to find ways to take into account the residual perturbative corrections and consider the question of whether the lattice results have perturbatively converged i.e., whether C n (z 2 ) is perturbatively convergent at all potential values of z that are used in a typical lattice analysis. There is an ongoing program to find higher loop corrections to C n , with the state-of-the-art result at next-to-next-to-leading order (NNLO) published in [45][46][47]. In a recent paper [48], the residual effect of higher-loops partially seen via resumming next-to-leading logarithms and threshold resummation was investigated in detail, wherein visible effects of resummation were seen in typical values of subfermi distances z. For the sake of argument, in addition to the resummation of lower-order logarithms, a priori it cannot be ruled out that the genuine higher-loop corrections and logarithmic terms that arise via 1PI diagrams could result in large numerical coefficients that multiply α m s factors, making the convergence of perturbation theory slower for the leading-twist approach. A method to take care of the renormalon ambiguity in the perturbation series was considered recently in [49].
In this paper, we assume two things: first, that the lattice data for the matrix element has a leading-twist description. Second, we assume that the Wilson coefficients that enter the leading-twist OPE are hadron independent. Given these two assumptions, the idea we pursue in this paper is to replace a perturbatively determined C n with a probabilistically likely C n that is determined directly from the lattice matrix element data, and given a trustworthy knowledge of valence PDF of a hadron from global analysis of experimental input. The lattice determined C n can then be used in computing the PDFs of other hadrons, as well as to obtain the zero-skewness GPD. By turning the argument around, as this Bayesian approach is falsifiable, it would also be a good way to rule out the hypothesis of higher-loop corrections being important in matching the lattice calculations (which is the premise of this paper), or whether there are substantial higher-twist corrections affecting the lattice results. Therefore, at the cost of adding the systematics from the lattice computations for two different hadrons, the method presented in this paper will help disentangle the perturbative QCD systematics from systematics of nonperturbative lattice parton computations. One should note that our method is a novel way to make use of global analysis data in lattice computation, unlike the previous ways [50][51][52][53], where one tries to include the lattice data as one of the input to global analysis. In the rest of the paper, we explain the method and apply it to lattice mock-data, and some real published proton and pion lattice data.

II. THE METHOD
In the traditional Bayesian reconstruction of PDFs f (x) via fits (of parameters that characterize a certain PDF ansatz) to the lattice data for the matrix element M H (ω, z 2 ) for a hadron H, one asks for the conditional probability distribution (4) One treats P (data|f ) ∼ e −χ 2 /2 via the χ 2 -minimization analysis to find the best fit value (i.e., the maximum a posteriori estimator) of f (x) that maximizes the lefthand side probability, and its confidence interval. The assumption of a certain parametrization of the PDF is implicit in the prior P (f ). Methods based on the Bayesian reconstruction of PDFs so as to avoid the modeling bias in f have been tried before [54]. Instead of applying the Bayesian methods to the reconstruction of PDFs, in this work, we apply the same ideas to the determination of Wilson coefficients. We propose the following two-step scheme assuming one has lattice results for more than one hadron -results for M H (ω, z 2 ) for a hadron H (say, for the proton), and the lattice data for M H (ω, z 2 ) for at least one additional hadron H (say, for the pion). A.
Step-1: Fit the Wilson coefficients given the lattice data for a hadron and prior knowledge of its pheno-PDF In the first step, we ask for the conditional probability distribution Q for the set of zand n-dependent Wilson-coefficients, themselves, given the lattice data for the matrix element M H (ω, z 2 ) for a hadron H, as well as the prior knowledge of the PDF f H (x) of the same hadron H from phenomenological determinations, As is standard, the distribution P (data for M H |{C n }; f H ) models the Gaussian fluctuation of the data around the theoretical model M twist−2 (ω, z; {C n }, f H ). That is, we can conveniently define, Q({C n }) = K exp −χ 2 /2 , with a normalizing constant K and where, σ H are the statistical errors associated with the lattice data for M H (ω, z). In our case, the theoretical model of the data is the twist-2 OPE that is dependent on the Wilson-coefficients as well the PDF (via its moments), where the sum runs over all z-values and all the values of order n in the twist-2 OPE that are included in the analysis. For example, one could use a broad prior width σ prior to input no bias, or another choice could be two or three times the variations coming from changing the scale from µ/2 to 2µ used to determine α s that enters C pert n (z). This first step to characterize Q({C n }) can be implemented in different ways depending on the computational effort that one is willing to invest. In perhaps the most accurate implementation, one can simulate the multidimensional integral, K e − 1 2 χ 2 ({Cn}) DC n , by using the Markov-Chain Monte Carlo (MCMC) numerical simulation. Thereby, one can estimate averages with respect to the distribution Q, of functions F ({C n }) that depend on the Wilson coefficients. The Expected A Posteriori (EAP) values of C n obtained by using F ({C n }) = {C n } is one such example that can computed using MCMC. A cheaper alternative to using MCMC is to first summarize the information in which can be obtained by the standard χ 2 -fits to the data by minimizing Eq. (6). The confidence intervals C MAP n ± σ Cn can be obtained by bootstrap. With the MAP values, a cheaper possibility to sample the distribution Q is to approximate it parametrically by the multinormal distribution N {C MAP n , σ Cn } . Instead of using a multinormal distribution, one could use the nonparametric bootstrap histogram of {C MAP n } obtained by performing the χ 2 minimization in each of the bootstrap samples of the data. For demonstration purposes, we will use this simpler procedure in the latter half of paper. In these approximations to the true Q({C n }), it is important to resample the central values C pert n entering the prior distribution from P ({C n }) in each bootstrap iteration, so as to ensure that in the absence of any constraint from the data, one reproduces the prior distribution P ({C n }) rather converging to the same value in every iteration. For the sake of completion on the possible implementations of Step-1, as a crude approximation, one could simply ignore the uncertainty in C MAP In the second step, we use the estimated C n from the first step to fit the M H (ω, z 2 )-data for another hadron H . That is, we ask for the conditional probability on the parameters α describing the PDF f H (x; α) of H given the data and marginalized over the random choices of {C n } picked from the distribution Q({C n }), The distribution, P (α), on the right-hand side is the prior on the parameters α that enter f H , given a fixed underlying model for the PDF. For example, f H (x; α) could be the model x a (1 − x) b , with α = {a, b}. In the rest of the paper, we will consider a non-informative prior on α, and hence, drop the prior term P (α) for the sake of simplicity. The distribution, P (data for H |f H , {C n }) ≡ K e −χ 2 (α,{Cn})/2 with the normalizing constant K , captures the Gaussian fluctuation of data for M H around the twist-2 OPE theoretical model via To be concrete, one is interested in a measure of the most probable value of f H and its confidence intervals. Similar to the discussion in Step-1, one such quantity that can be measured in an MCMC simulation is the EAP value of the PDF, To evaluate the above integral, one will sample the values of {C n } from the MCMC simulation as performed in Step-1 in order to capture the measure Q({C n })DC n , and perform a second MCMC simulation over the parameters α. This would be computationally rather complex. Instead, as a simpler statistic for the probable value of , and then perform an average over C n . That is, one can use as an estimator of the PDF and its confidence intervals.
To perform the average . . . Cn , one needs to draw instances of C n sampled in Step-1. As mentioned previously, we will use the bootstrap sample of {C MAP n } from Step-1 as an approximation of Q({C n }) to perform the above step. To avoid extra symbols, we will simply use f H (x) to denote f avg H (x) and its error-bands in the rest of the paper.
The two steps could be done in such a way that collaboration-A does a dedicated computation of an ensemble of {C n } distributed according to Q({C n }) in its own set of gauge ensembles etc., using some hadron H, and collaboration-B takes this Q({C n }) as an input for its dedicated computation of f H (x) for hadron H . In this case, there is no correlation between a choice of C n and the re-sampled bootstrap configurations. In another way, a single lattice collaboration computes both pion and proton, preferably at the physical point and using the same set of configurations. In this case, one could reduce the statistical errors by sampling Q({C n }) so as to maintain correlations between the pion and proton data.

C. Remarks on the method
First, we have a choice as to which hadron we should consider first to determine the Wilson coefficients. If one has very accurate data, this question simply boils down to how well determined the pheno-PDFs of a given hadron is. For the unpolarized sector, using the proton valence PDF, which has been determined at NLO and NNLO level accuracy [39] from the experimental data, as the known input would be the ideal choice. One can then input the set of C n from the proton analysis into the analysis of pion PDFs on the lattice. The downside to this order of steps is that the moments x n of proton PDF approach zero faster as n is increased, compared to the pion 3 . Consequently, while one can determine C n (z) for to a higher order in n (we expect about up to first three even n, as we will discuss later in the paper) using pion PDF as input, it would be difficult to do so using the proton. Nevertheless, this difficulty is purely statistics related and can be solved in a high-statistic study of the proton ITD in the future. On the other hand, the difficulty in using valence pion PDF as the input PDF is that unlike the proton PDF, the robustness of valence pion PDF from different analyses of experimental data [40,[55][56][57][58][59][60][61][62][63][64] are still being scrutinized. In this paper, to simply demonstrate the idea, we will use the pion PDF as the starting point, but use different versions of pheno-PDFs of pion as inputs in the lattice analysis, and then determine the C n for each such choice. In turn, this procedure will lead to different estimates of proton PDFs on the lattice. In this way, in addition to demonstrating the method, we would also be studying the impact of a good phenomenological determination of pion-PDF on the lattice determination of proton PDFs.
Secondly, let us elaborate on the method being a Bayesian extension of the traditional analysis of lattice PDFs. As our confidence in the perturbative Wilson coefficients increases, due to say its calculation up to a certain higher-order m, then by the Bayesian framework, one will choose the prior distribution that peaks more and more about the m-loop result for {C n }. That is, one will choose σ prior → 0 as the confidence is subsequently increased. In that case, the posterior distribution approaches a delta-function about the m-loop result, and one recovers back the usual way of analysis using a given perturbative C n . This is the essence of the idea in this paper to treat our ignorance of C n in a probabilistic manner.
Thirdly, it is easy to see that the fitted coefficients will automatically satisfy the Altarelli-Parisi evolution to the same accuracy to which the input pheno-PDF moments satisfies the evolution, since the matrix element M H , that is renormalized in a lattice scheme such as the ratio, is independent of the factorization scale µ. To expand on, the actual quantities that we obtain by fits are the z-dependent coefficients of the (−iω) n /n!terms in the twist-2 OPE in Eq. (7). Let us call those coefficients D n (z). Since M H has no knowledge of µ, the coefficients D n (z) will be constants independent of µ. By writing D n (z) = C n (µ, z) x n (µ), it is clear that the µ-dependence of fitted C n (µ, z) will be the inverse µdependence of x n (µ), which is dictated by the Altarelli-Parisi equation that the global fit PDF satisfies.
Fourthly, since we are using only the valence Mellin moments in the twist-2 OPE to determine C n , it is not necessary to use the pheno-PDF values, and instead, one could use lattice determinations of Mellin moments obtained using the local twist-2 operators (e.g., [65,66] for recent studies). Using the local twist-2 operator approach, both even and odd Mellin moments can be computed. However, currently, only x 2 even valence moment from such studies are usable as input to determine C n in the method we proposed (and only even n enter the valence ITDs), and hence not a practical option currently.

III. GENERALIZATION TO ZERO-SKEWNESS GPD
Let us now consider the equal-time matrix element corresponding to the generalized parton distribution function (GPD) f (x, t) at zero-skewness [67,68], with spatial momentum P = (0, 0, P z ) and transverse momentum transfer q = (q x , q y , 0), with t = q 2 . For this kinematics at zero skewness, ξ = q z /2P z = 0, the corresponding leading twist OPE is given as [67,68] , with the same set of Wilson coefficients {C n } as in the case of forward PDF [67] 4 , but the GPD moments (i.e., the generalized form factors), x n (t, µ) being dependent on the momentum transfer t. Given this property of the ξ = 0 GPD, it is straight forward to generalize the analysis scheme presented in the last section: 1. Determine Q({C n }) using the lattice forward matrix element data M H (ω, z 2 , t = 0) of a hadron H using the pheno-PDF f H (x) as a given. This step is the same as the step-1 described previously. The data from z ∈ [a, 12a] have been put together in the plot (and not distinguished by their z values to avoid clutter), leading to visible corrections to the scaling with ω and will be captured by the Wilson coefficients Cn(z). For parameters used to generate this data, refer text.
We expect this to be a less-noisy implementation as the GPD data of a hadron H is expected to be correlated with the t = 0 matrix element of the same hadron H in a given ensemble. It would be interesting to see the data presented for valence pion GPD [69] and unpolarized proton GPD [70] analyzed in the manner suggested here.
In addition to the generalization to zero-skewness GPD, few other cases are also possible 5 . The chiralsymmetry present in the massless quark limit ensures that the Wilson coefficients are the same for bilocal quark bilinear operators with γ z and γ 5 γ z Dirac spinor structures [35,71]. That is, the matching factors are the same for the unpolarized PDF computed with γ z bilocal quark bilinear and the helicity PDF computed with γ 5 γ z bilocal bilinear. Thus, the method presented here can be applied to the proton unpolarized PDF determined using γ z bilocal bilinear to determine the corresponding {C n }, and then use them in the lattice computation of the proton helicity PDF. The drawback is that the γ z bilocal bilinear suffers from mixing with the scalar bilinear [43,72], adding to the systematics. In addition, the Wilson-coefficients for γ z bilocal bilinear also enters the leading-twist OPE that describes the Distribution Amplitude (DA) of the pion [9,68,73].

IV. DEMONSTRATION OF THE METHOD USING MOCK-DATA
In this section, we apply the idea to a set of mock-data that we generated for the pion and proton pseudo-ITD M(ω, z 2 ), which are renormalized in the ratio scheme and 5 We thank the anonymous referee for pointing us to the additional possibilities.
set in a lattice with a spacing a = 0.06 fm with an extent L z = 48 in the z-direction. Also, we will only be looking at the real part of the ITD, ReM, which corresponds to the nonsinglet valence contribution for both the proton and the pion. The point of using mock-data is that we exactly know what corrections to NLO matching we put in, and ask how this new framework will handle these corrections. We generated the mock-data for pseudo-ITD of H (proton p and pion π) by drawing sample means M H (zP z , z), for z ∈ [a, 12a] and P z aL z /(2π) ∈ [1,5].
where the moments x n H were inferred from the pheno-PDF values of H; namely, from the unpolarized NNPDF [39] data at µ = 3.2 GeV for the proton, and the JAM20 [40,64] data at µ = 3.2 GeV for the pion respectively. Since we are looking at the real part of the ITD, only the even n-values enter the above equation. We use a set of mock Wilson coefficients C mock n given by the coefficients of the ln j (µ 2 z 2 ) terms, c n (j) =c NLO n (j) + ∆c n (j), j = 0, 1, 2, used in Eq. (3). Above,c NLO j are the NLO values that can be read off from [27,44], and ∆c j are toy correction terms that we put in by hand to simulate the higher-loop corrections. Since we are working with ratio matrix elements, we simply take C mock 0 (µ 2 z 2 ) = 1, as it is true to all orders. We added the residual higher-twist correction via the term in the second line of Eq. (17); it is motivated by its expansion as Λ 2 z 2 ω 2 + Λ 4 z 4 ω 4 /2 + . . ., which involves all z k type terms, and we have chosen it as a function of z 2 ω 2 as any leading z 2 ω 0 type highertwist effect would most-likely have been canceled in the ratio scheme. As for the covariance matrix Σ in the multinormal distribution, we made an ad hoc choice Σ [z,ω],[z ,ω ] = σ z,ω σ z ,ω /(1 + k|ω − ω | δ ), with the diagonal term being simply the σ 2 z,ω error term, and the |ω − ω | δ term controlling how the correlation in the mock-data falls off as one moves along the off-diagonal elements; we chose δ = 0.5 and k = 0.5 for maintaining some correlation in the mock-data. As for the diagonal error term, we chose a model σ z,ω = r|ωz| motivated by the fact that for a given value of P z , the data gets noisier as z increases, and that at a fixed z, the data at a higher P z is noisier. We used r = 0.0005 for the pion and 0.0008 for the proton. The only justification for choosing all these models to generate mock-data is simply that it leads to reasonable looking data that is similar to what is seen in the published literature.

A. Case-1: When the dominant corrections are higher-loop corrections to leading-twist OPE
First, let us look at a case where there are only corrections to the NLO C n , which we will denote as C NLO n , via non-zero ∆c n (j) terms, and no higher-twist correction terms, i.e., Λ = 0 for both the pion and the proton. This is exactly the case where the our method has to work better than the one using finite order C n , and therefore, it should not be surprising that we will account for the corrections at the end of this subsection. We experimented with different models of corrections to NLO Wilson coefficient (and the reader can try their own setup in the Jupyter notebook here [74]), but we discuss an example case that we think is more plausible, with corrections given as [∆c n (0), ∆c n (1), ∆c n (2)] = [0.1037, 0.0207, −0.0025] ln(n + 1).
The choice of ln(n + 1)-dependence of the correction terms is inspired by the asymptotic behavior of the Harmonic functions that appear in the perturbative kernel and could be important to threshold resummation as discussed in [48]. The mock-data for ReM that is generated from N M 0 H , Σ for the proton and the pion are shown in Fig. 1.
By using the JAM20 [40,64] determination of pion PDF as the fixed input, we fitted the coefficients C n (z) for even values of n > 0 as the free parameters to the mock-pion data using the twist-2 expression in Eq. (7) with the truncation at n = N max = 10. As discussed above, for the ratio pseudo-ITD, the value of C 0 (z) = 1, and hence it is not a fit parameter. The fits are analogous and non-parametric in z to the type of fits used in "OPE without OPE" analysis [44] of moments as a function of z, with the modification being that we are now fitting the z-dependent coefficients C n (z) instead of moments. We imposed a prior on C n from C NLO n with a rather broad width of σ prior of 100% of C NLO n at any z. In the top panels A, B and C of Fig. 2, we show the three of the fitted C n (z) for n = 2, 4 and 6 as a function of z. These central (MAP) values and their errors shown in the plot (and the covariances between different C n (z) which are not shown) summarize the abstract posterior distribution Q({C n }) we used in our discussion of the method. For comparison, C NLO n is also shown as the black curve, which clearly disagrees with the fitted C n , especially for C 2 . This is expected by construction in this example mockdata. For z/a = 1 for C 4 and z/a ≤ 4 for C 6 , the fitted values fall exactly on the NLO curve. This is due to the prior imposed, as there is no information on C n for these z and its knowledge is purely derived from the prior. As z is increased, one sees a cross-over from a prior-driven fit to a more data-driven one.
In the second step, we used the C n obtained above in the twist-2 OPE analysis of the mock proton pseudo-ITD data shown in Fig. 1. We performed a variant of the "OPE without OPE" analysis, in which we fitted proton moments x 2 p , x 4 p and x 6 p as free fit parameters in the twist-2 OPE expression, to the proton mock-data satisfying, z ∈ [2a, z max ] and all momenta. We took care of the uncertainty in C n by randomly drawing from the bootstrap samples obtained in the analysis of the pion in the first step. We repeated these fits by increasing z max gradually. We have shown the dependence of the fitted moments x 2 p and x 4 p on z max in the bottom panels D and E of Fig. 2, respectively. As expected, x 2 p from the fits using C n obtained from the pion agrees with the NNPDF value, which is the underlying PDF used to generate the proton pseudo-ITD mock-data. Similar to the case in actual lattice computations presently, the signal for x 4 p is poor at shorter z max , and albeit up to errors, the results with fitted C n agrees better with NNPDF value than the NLO C n as we venture to larger z max . In the bottom panel F of Fig. 2, we demonstrate how this method can also be applied in the reconstruction of xdependent PDF by using a two-parameter ansatz of the form, f (x) ∼ x α (1−x) β ; for this fit, one simply writes the moments as a function of α and β, and then performs the fits using the twist-2 OPE formula in Eq. (7) with {C n } being from either NLO or from the fits to the pion. The PDF reconstruction analysis with fitted C n reproduces the proton PDF, f p (x) to a better accuracy, whereas the one using the NLO kernel does not, as expected in this example. Through the above analysis of the mock-data, we have elaborated the typical steps we expect to be involved in the implementation of the method.

B. Case-2: Possibility of suppressing higher-twist effects
The method we have laid out still assumes that the leading twist OPE describes the lattice data to a good accuracy. However, given that we only know the data up to errors, we can ask if the higher twist corrections can effectively be absorbed into effective Wilson coefficients determined by fits; we have to admit that this procedure is not pristine. In order to look for this possibility, let us consider a higher-twist correction to Eq. (7) of the form Λ 2 π z 2 ω 2 for the pion pseudo-ITD and Λ 2 p z 2 ω 2 for the proton pseudo-ITD. Let us first assume that such a higher-twist effect is more or less hadron mass independent, and that Λ π ≈ Λ p = Λ. In the first step, when we do the fits of {C n } to, say the pion, such a higher-twist term can effectively be absorbed into a Wilson-coefficient given by where C twist−2 2 (µ 2 z 2 ) is the part of the fitted value that is the actual twist-2 Wilson coefficient. When such an effective C eff 2 (z) is used in the analysis of the proton pseudo-ITD, the net higher-twist correction in proton gets re-  (Bottom panels D, E) The Mellin moments, x 2 and x 4 , extracted by twist-2 OPE fits to the proton pseudo-ITD over the range of z ∈ [2a, zmax] are shown as a function of range maximum zmax in the panels D and E respectively. The results using NLO Cn(z 2 ) as well as the Cn(z 2 ) obtained from fits to the pion ITD (shown in the top panels) are displayed together. The corresponding NNPDF moments that were used to produce mock-proton data are the black horizontal lines. In bottom panel F, the results from fits of the JAM-type two parameter ansatz f (x) = x α (1 − x) β are shown using the two types of Cn.
duced to Λ 2 1 − x 2 p x 2 π z 2 ω 2 , which by putting in values for the moments at µ = 3.2 GeV, can be estimates to be a 70% reduction. On the other hand, if the higher-twist correction term Λ p , Λ π were proportional to the proton and pion masses respectively, there would hardly be any reduction in higher-twist term for the proton by this process. Unfortunately, the property of the higher-twist correction in the lattice data are unknown.
We repeated the mock-data analysis presented in the previous sub-section, but this time, included the same value of Λ = 0 in Eq. (17) for both the pion and the proton. We used a value of Λ ≈ 30 MeV, which was found in the pseudo-ITD in RI-MOM scheme [30], and even smaller in the ratio scheme. We performed the analysis in the similar way as described in the previous subsection by inputting the JAM20 pion-PDF to determine C n (which are effective now due to non-zero Λ). We show the result of the reconstructed proton PDF in Fig. 3. The use of fitted C n brings the estimated PDF closer to the actual proton PDF underlying the data when compared to the result obtained by using the NLO C n . But, the results still disagree with the NNPDF curve due to the fact that the higher-twist corrections are canceled only partially.
x Mock data with HT NLO C n fit C n NNPDF   FIG. 3. The result of x-dependent PDF fit to the proton mock-data containing higher-twist (HT) correction (see text) using a two parameter ansatz f (x) = x α (1 − x) β . The results using NLO Cn and using an effective Cn from fits to the mockdata for the pion are shown. The actual underlying PDF (NNPDF) that is used to produce the mock-proton data is shown by the black curve.

V. APPLICATION TO SOME SELECTED PUBLISHED RESULTS
In the last part of the paper, we apply our methodology to some of the previously published lattice results from different collaborations. However, we should remark  [30]. The fits were made using two different pheno-PDFs at µ = 3.2 GeV, namely, using the JAM20 PDF (red) and the ASV soft-gluon resummed estimate (blue). The NLO results for Cn are the black curves. Bottom panels D, E, F show the application of Cn extracted from the pion, to the analysis of the proton pseudo-ITD data presented in M. Bhat et al [34]. The panels D, E show the fitted values of x 2 p and x 4 p as function of the maximum value of z = zmax of the fit range z ∈ [2a, zmax]. The red, blue and the black points correspond to the estimates using the values of Cn obtained using JAM20 pion PDF, using ASV PDF and the NLO values respectively. The estimate of the moments from the NNPDF at µ = 3.2 GeV are shown as horizontal black bands. The panel F show the proton PDFs reconstructed using two-parameter ansatz using the three-types of Cn. that our analysis is certainly crude as we simply take the central values and errors presented in the publications, and we do not take into account the correlated statistical fluctuations in the data. More importantly, we will apply this method to some cases where the computations were performed with heavier than physical pion masses. We believe this is not a serious issue, at least for the pion PDF (c.f., [26,75]). However, looking forward to the future, where more computations of both the pion and proton at the physical point might be forthcoming, we expect this method will be more appropriate. However, with these cautions, let us apply the method to some of the existing data in the literature.
The details of the published lattice QCD results for the pion and proton that we borrow the data (centralvalue and error) from, are as follows. For the pion, we take the ratio pseudo-ITD results [30] for a 300 MeV pion determined on an ensemble with a = 0.06 fm and lattice size of 48 3 × 64. For the proton, we take the ratio ITD results at the physical point [34] obtained on an ensemble with a = 0.0938 fm on 48 3 × 96 lattices. For complete details on the lattice methods used in these studies, we refer the reader to the respective papers [30,34].
We followed the same set of steps as described for the analysis of mock-data sets for these real lattice data, and hence we will keep the details brief. We used the pion PDF as the starting point. However, in the analysis here, in addition to the JAM20 pion PDF, we also included the ASV pion PDF estimate [63] (obtained by a reanalysis of the Fermilab data [57] by including soft-gluon resummation). In the pion PDF literature, there is still an ongoing debate on whether the large-x exponent β is closer to 1 or about 2 [63,[77][78][79][80][81][82][83]. In the top panels (A, B and C) of Fig. 4, we show the Wilson coefficients C 2 , C 4 and C 6 obtained by fits to the pion pseudo-ITD data. The fitted values of C n obtained by assuming the JAM20 PDF as the input PDF, are shown as red circles in the plots. On the other hand, when we assumed that the ASV pion PDF at the same µ = 3.2 GeV scale as the given input PDF, we obtained the points shown in blue. For brevity, we will refer to the two sets of Wilson coefficients as C JAM n and C ASV n , respectively. We remark that it would be a good idea to skip the first one or two lattice points for C 2 in the above plots as they are likely to be affected by small-z lattice artifacts. To compare, the NLO values of C n are shown as the black curves, which lie comparatively closer to C JAM n than to C ASV n . To rephrase, if one is confident that the ASV result is the correct pion PDF, then it implies that even the first non-trivial Wilson coefficient C 2 needs to have about 40% corrections to NLO obtained from fit to pion pseudo-ITD data (taken from X. Gao et.al [30] and shown in the top panels of Fig. 4) are applied to a completely different pion pseudo-ITD result taken from B. Joó et al. [76] to check for internal consistency and weak dependence on mass and lattice spacing. The result of two sets (red and blue) of fitted moments x 2 π obtained from the usage of C JAM  [33] in the same ensemble as [76], and shown as a function of zmax. The NNPDF result is shown for comparison. C 2 resulting from higher-loop terms (given the assumption that higher-twist terms are not as important). For C 2 and C 4 , the fits are driven by the lattice data as well as the input PDF for all z, whereas for C 6 , the data is driven only by NLO prior for z < 0.3 fm and then there is cross-over to data-driven regime for larger z.
In the next step, we applied C JAM n and C ASV n determined above to the proton pseudo-ITD lattice data. We dealt with the mismatch of lattice spacing between the two lattice calculations by interpolation. We first fitted the even moments, x 2 p and x 4 p to the proton pseudo-ITD data using twist-2 OPE truncated at O(ω 4 ) containing the choices of Wilson coefficients as C JAM n , C ASV n and C NLO n at µ = 3.2 GeV. We performed these fits including the data with z ∈ [2a, z max ] for z max up to 0.75 fm. The results for x 2 p and x 4 p from these fits are shown as a function of the maximum of the fit range, z max , in the panels D and E of Fig. 4. The black, blue and red points are the results from such fits using C NLO n , C ASV n , C JAM n respectively. For comparison, the black bands are NNPDF estimates of the proton moments at µ = 3.2 GeV. For z max ≈ 0.24 fm, there is a discernible signal for x 2 p , while x 4 p remains comparatively noisier. To account for the possibility of an underestimation of the errors in our analysis due to the absence of any correlated statistical fluctuations, we conservatively included both the 1-σ as well as the 2-σ errors in the plots for the moments. The NLO result for x 2 p is a bit higher than the NNPDF value (as also observed in the actual full-fledged analysis in [34]). The usage of C JAM n leads to a slight decrease in the estimated value of x 2 p , and surprisingly, the largest change comes from the usage of C ASV n , which lowers the estimated x 2 p that also agrees better with the NNPDF value. For x 4 p , there are corresponding changes in the central values depending which C n is used, but completely masked by the larger error-bars. In the bottom panel-F, we show the unpolarized proton PDF, f p (x), that we reconstructed using a simple, f p (x) ∼ x α (1 − x) β ansatz. Again, the results obtained using the three types of C n are shown in the different colors; the 2-σ estimate is the outer-band and 1-σ result is the inner-band enclosed by the solid lines. The x-dependent results are sensitive to which Wilson coefficient is used, and particularly, the result obtained using C ASV n agrees better with the NNPDF curve. In order to see if the results presented above are some artifacts of putting together the results from two completely different lattice calculations, at different pion masses (with only the proton at physical point) and at different lattice spacings, we repeated the above analysis to the central-value and error data of the pion and proton pseudo-ITDs that we borrowed from the papers [33,76] by the HadStruc collaboration. The results reported in these papers used 32 3 × 96 lattice with a = 0.127 fm at a heavier pion mass of 415 MeV. The analysis of the 415 MeV pion data from Ref. [76] by applying the set of C n obtained also from the pion, with a 300 MeV mass and from a completely different data in [30], is an interesting exercise. If there is not much lattice spacing and pion mass dependence in the two pion data, then the analysis using C JAM n from panels-(A,B,C) in Fig. 4 should result in the JAM20 moments as it is the initial input PDF in this case, and similarly, an analysis using C ASV n should result in ASV moments. In the panel-A of Fig. 5, we show x 2 π from such an analysis by fitting the moments over a range z ∈ [2a, z max ]. We have again conservatively plotted both the 1-σ and the 2-σ errors. It is satisfying that this pion analysis reproduces the JAM20 and ASV moments as per our initial expectation. In the right panel-B of Fig. 5, we similarly show the result of applying C JAM n , C ASV n and C NLO n to the proton data from [33]. The result is very similar to the one we observed in the panels-(D,E) of Fig. 4, with a similar effect of the choice of C n on the determined moments.
Given the cautions we offered at the beginning of this section, we need to be careful about the conclusions we can draw from the observations we have made. First, one should pay attention to fact that the usage of ASV versus JAM20 pion PDF as input PDFs leads to changes in C n (z) even for n as small as 2. To contrast, in Ref. [48], the corresponding threshold resummation for the pseudo-ITD calculation was investigated as the effect of potentially large α s ln(n) and α s ln 2 (n) terms that will be present in C n at larger n, and thus will not affect the smaller n. Thus it is paradoxical that imposing the correctness of ASV PDF automatically implies a big correction to C n even for small n; we do not have a resolution in this paper. Thus, it would be important to repeat the ex-ercise presented in this section to the lattice data for pion pseudo-ITD at the physical point when available in the future to see if the observation we have made still stands. Secondly, we see that an unaccounted larger higher-loop correction to C n starting from n = 2 is sufficient for a better agreement between lattice computations of proton PDF and the phenomenological estimates, without any necessity to invoke higher-twist contaminations to the leading-twist OPE. The curious agreement of the lattice estimation for the proton moments and the PDF when C ASV n is used in the analysis instead of NLO C n , brings forward the great impact a better determination of the pion PDF from global analysis in the future (c.f., [84,85]) will have on a lattice computations of proton PDFs; a somewhat unexpected and interesting connection due to the method presented here.

VI. DISCUSSION
The method we advocated in this paper for lattice determination of valence PDFs consists of the following steps: (1) The leading-twist Wilson coefficients C n (µz) are determined directly by fitting the lattice data for a hadron, say the pion, by assuming the pion PDF is known well from phenomenological determinations and then taken as an input for the lattice analysis. At this step, our degree of confidence that a perturbative calculation at a given order is sufficient can be incorporated via addition of priors to the fits. (2) The C n (z) obtained from fits is then used in the analysis of lattice data to determine the PDF of another hadron, say, the proton by performing the usual reconstruction analysis methods to determine the x-dependent PDF, except that one replaces the perturbative matching coefficients with the one obtained from the pion, in this example. (3) One can also use the fitted values of C n in the analysis of zero-skewness GPD of, say the pion itself, at non-zero momentum transfers. We demonstrated the feasibility of this method by applying it to a set of mock-data, as well as to some published lattice pseudo-ITD data for the pion (at 300 MeV and 415 MeV) and the proton (at physical point and at 415 MeV). In doing so, we showed how the precise extraction of pion-PDF from the experimental data can have an impact on lattice studies of proton PDF, indirectly via the method we presented.
A drawback of the method is that there should be a good determination of pheno-PDF; this is achievable for the valence unpolarized PDF, but fails for the nonsinglet helicity and transversity PDFs since they are defined only for the proton. In such cases, the analysis strategy (3) is still viable, i.e., extend to zero-skewness helicity and transversity GPD of proton by inputting the experimental values of helicity and transversity PDF of proton respectively. This can also be tested in the helicity GPD data already presented in [70]. In the future, it would also be interesting to extend this study to the iso-singlet PDF cases wherein the small-x uncertainties in the global fit PDF themselves might need to be propagated into the analysis of lattice matrix elements (for example, see Ref. [86] for a discussion of uncertainties from small-x logarithms in global fit determination of singlet PDFs).
One could ask why should one perform this rather complicated procedure when one is, say, 75% certain that the correction to the Wilson coefficient C 2 (z) from higherloop calculation will be less than 2%. In fact, the point of the paper was to generalize the perturbative matching method in such a way as to fold in our uncertainties about perturbative Wilson coefficients itself in the analysis; the Bayesian approach achieves this for us. In the example, we would arrange the Bayesian prior distribution of C 2 (z) in such a way that it is 75% likely that the C 2 lies within 2% of the NLO Wilson coefficient. This way, we are making the degree of our ignorance of C n quantitative.
A reasonable concern could be that the experimental data is known only up to NLO or NNLO order, and hence, why should one use a procedure other than using NLO or NNLO value of C n for matching. The assumption we are implicitly making here is that the perturbative convergence of pheno-PDFs is fast and wellunderstood, whereas it does not necessarily guarantee a rapidly converging perturbative Wilson coefficients that are used to match PDFs to the lattice matrix elements. Thus, the solution, which in retrospect appears very simple, is to replace the two different perturbative calculations (one to extract pheno-PDF and the other to match lattice results), with a single well-studied perturbative step that has been used to determine the pheno-PDF via global analysis.