Extraction of Next-to-Next-to-Leading-Order PDFs from Lattice QCD Calculations

Quark correlation functions in position space are calculable in lattice QCD and factorizable into parton distribution functions with matching coefficients perturbatively calculable to all orders in QCD, which provides a way to extract PDFs from lattice QCD calculation. We present for the first time complete next-to-next-to-leading-order calculation of valence-quark matching coefficients. We find that theoretical uncertainties are improving with higher order contributions. Our method of calculations can be readily generalized to evaluate sea-quark matching coefficients and gluon correlation functions, putting the program to extract partonic structure of hadrons from lattice QCD calculations to be comparable with that from experimental measurements.


INTRODUCTION
Parton distribution functions (PDFs) encode important nonperturbative information of strong interactions, and they are crucial for understanding all phenomena at the Large Hadron Colliders (LHC) [1]. In terms of QCD factorization [2], a typical hadronic cross section with a large momentum transfer Q and collision energy √ S at the LHC can be factorized as where i, j = q,q, g represents parton flavor, f i/h (x, µ 2 ) is the PDF as a probability distribution to find an active parton of flavor i inside a colliding hadron h with the parton carrying the hadron's momentum fraction x, probed at a factorization scale µ ∼ O(Q), dσ ij represents a short-distance partonic scattering, and ⊗ indicates an integration over value of x or x , accessible by the scattering cross section. By measuring hadronic cross sections, with perturbatively calculated partonic hard parts dσ ij , PDFs have been extracted from the world data at the state-of-the-art next-to-next-to-leading order (NNLO) accuracy [1]. With the steep falling nature of PDFs as x → 1 and the convolution in Eq. (1), the uncertainty of extracted PDFs at large x is so significant that limits our confidence to push the search for signals of new physics to larger invariant mass. With the nonperturbative nature of PDFs, it is natural to ask if we can calculate PDFs directly in lattice QCD (LQCD). A short answer is no since the operators defining PDFs are time-dependent and LQCD is formulated in Euclidean space-time. Recently, stimulated by the seminal work of quasi-PDFs approach proposed by Ji [3], extraction of PDFs from lattice QCD calculation has drawn many new ideas, including the pseudo-PDFs [4], current-current correlators in momentum space [5] and current-current correlators in position space [6]. See also some earlier related approaches [7][8][9][10][11][12].
As proposed by two of us in Refs. [6,13], PDFs can be extracted from any good LQCD observables, which was referred to as "Lattice Cross Sections" (LCSs), that are calculable in LQCD and factorizable into PDFs with perturbatively calculable matching coefficients, where ξ with ξ 2 = 0 represents the size of nonlocal operator O n (ξ) of type n, controlling the short-distance physics of the factorization, and ω ≡ p · ξ (often referred as Ioffe time). With perturbative matching coefficients K n/i between σ n/h (ω, ξ 2 ) and f i/h (x, µ 2 ), PDFs can be extracted by QCD global fits of data generated by LQCD calculation of σ n/h (ω, ξ 2 ) with various operator type n, just like how PDFs have been extracted from the world data on various high energy scattering cross sections [1,6,13]. σ n/h in Eq. (2) can be parton correlation functions defining quasi-PDFs [3], current-current correlators [6], or any other quantities that satisfy the aforementioned properties. Since σ n/h does not have to be a physical cross section, the operator O n (ξ) defining σ n/h might require additional ultraviolet (UV) renormalization beyond using renormalized fields. This additional UV renormalization has impacted perturbative calculation and stability of the matching coefficients K n/i for LQCD observables. Although extraction of PDFs from LQCD calculations have made tremendous progresses in recent years , the state-of-the-art calculation of short-distance matching coefficients is still limited to the next-to-leading order (NLO) in almost all existing approaches [27][28][29][30][31][32][33][34], which is partially limited by this additional renormalization and our ability to do perturbative calculation in co-ordinate space. In this Letter, we derive for the first time the NNLO valence-quark matching coefficients in dimensional regularization, allowing us to extract PDFs from LQCD calculations at the same rigor as those extracted from experimental data. Calculating the matching coefficients at NNLO is not only necessary for being competitive with the effort to extract PDFs from experimental data, but also for addressing concerns that a possible mechanism may break the factorization at NNLO [54].

QUARK CORRELATION FUNCTIONS
We focus on the following gauge invariant quark correlation operator [3] which is made of renormalized fields with a path ordered gauge link in the fundamental representation, Because this composite quark correlation operator is UV divergent, a UV regulator δ is needed, which may represent lattice spacing a in lattice QCD calculations, or represent ≡ (4 − d)/2 in dimensional regularization (DR) of continuum calculations. µ is a dimensional scale accompanied by the UV regulator, which is different from the factorization scale in Eq. (2), while one could choose them to be equal numerically. This UV divergence is multiplicatively renormalizable [21][22][23], as where superscript RS indicates a renormalization scheme and Z RS (ξ 2 , µ 2 , δ) is the multiplicative renormalization constant. We will choose regularization-invariant renormalization conditions so that the renormalized O ν,RS q are independent of δ and µ 2 .
Quark correlation functions (QCFs) are defined as hadronic matrix elements of O ν,RS q (ξ) which is independent of regularization scheme and scale, like physical cross sections. With ξ 0 = 0 and ξ 2 Λ 2 QCD 1, F ν,RS q/h (ω, ξ 2 ) are expected to be calculable in LQCD and factorizable into PDFs [3] with an all-order proof of the factorization given in Ref. [13]. Their Fourier transform over dω with fixed p leads to the quasi-PDFs; and with fixed ξ is proportional to pseudo-PDFs [6]. In this Letter, we focus on flavor non-singlet case, and have corresponding factorization formula in continuum as [6] is a finite renormalization factor that transforms any "preferred" regularization-invariant RS scheme to the conventional MS scheme, K ν are perturbative matching coefficients in MS scheme, and q v ≡ q −q means To extract the valence quark distribution f qv/h from LQCD calculations of F ν,RS qv/h to the NNLO accuracy, we have to perturbatively calculate R RS and K ν to the power of α 2 s .

RENORMALIZATION CONSTANT
The renormalization constant Z RS introduced in Eq. (4) is determined by short-distance property of the quark correlation operator in Eq. (3) and should not depend on the hadronic state used to define the QCFs of this operator. Because of its multiplicative renormalizability, matrix element of O ν,b q in Eq. (3) with any state could define an allowable renormalization scheme, wheren is any vector keeping the denominator nonvanishing and the superscript "(0)" indicates that the matrix element is evaluated to the lowest order in perturbation theory. Different choice of the state |RS corresponds to different renormalization scheme. For example, an offshell quark state with a specific momentum was used in defining RI/MOM or RI /MOM scheme [24][25][26][27][28]; a hadron state with zero momentum was used in calculations of pseudo-PDFs [33] [Matrix element in this case cannot be perturbatively calculated and one may choose the denominator in Eq. (9) as 1]; and the vacuum state was introduced in Ref. [55].
In the following, we define the renormalization constant with the vacuum state |Ω and denote RS = vac. By calculating the vacuum expectation value to NNLO, we demonstrate that without an identified external momentum, the renormalization constant Z vac is completely free of infrared (IR) and collinear (CO) singularity and its UV divergence is regularized by DR, from which we obtain Z MS (ξ 2 , µ 2 , ) and R vac (ξ 2 , µ 2 ) at NNLO level.
In Fig. 1(a,b,c), we show some representative Feynman diagrams, up to NNLO, for the vacuum expectation Ω|n· O b q (ξ, µ 2 , δ)|Ω . The diagram (a) in Fig. 1 determines the normalization of Z vac , where |ξ| 2 ≡ −ξ 2 , and the result agrees with Ref. [55]. The Fig. 1 where we assume without loss of generality that zcomponent ξ z is only nonzero component of ξ, andn satisfiesn · l ≡ l z for any vector l. We find that it is convenient to carry out the integration in Eq. (11) by Fourier transforming the ξ z into q z in momentum space as F [M b ] ≡ dξ z e −iξzqz M b to eliminate the exponential factor by using where 2πδ(x) = −2 Im( 1 x+i0 + ) is used. The Fourier transformation also ensures that only imaginary part of gauge-link-related propagators are involved, which led to a similar effect of optical theorem. Our matrix element is defined with gauge-link in coordinate space, which is effectively equal to sum over diagrams with cut gaugelink in momentum space. It is the summation of cuts of gauge link that forces the appearance of imaginary part of "forward scattering amplitude". The obtained loop integrals in momentum space can be reduced to linear combination of a small set of integrals, called master integrals (MIs), by using integration-by-parts relations (IBPs) [56,57]. We use the package FIRE5 [58] to do this reduction, which results in with two vacuum MIs defined as , To carry out these single-scale vacuum MIs, we use the method presented in Ref. [59] by setting up and solving dimensional recurrence relations and obtain We then Fourier transform inversely from q z dependence into ξ z dependence to derive the result of M b in DR.
Other two-loop diagrams, including UV counter term diagrams, can be calculated similarly. All three-loop diagrams like diagram (c) in Fig. 1 can also be calculated similarly as the diagram (b) described above. The only difference is that analytical expression of vacuum MIs cannot be obtained by solving dimensional recurrence relations directly. Instead, we calculate the vacuum MIs to high accuracy by using dimensional recurrence relations and then obtain exact results by using PSLQ algorithm [60]. We check the correctness of our exact results numerically with at least 10 3 digits.
By adding all diagrams and UV counter terms together, the remained divergences should be removed by operator renormalization. With a MS subtraction scheme, we obtain Z MS (ξ 2 , µ 2 , ) and R vac (ξ 2 , µ 2 ) at NNLO level, with analytical expressions given in Appendix.

MATCHING COEFFICIENTS
By choosing the MS scheme for QCFs, we have the same factorization in Eq. (6) with R RS = 1 and introducing a µ dependence on the left hand of the equation. To calculate the matching coefficients K ν , we replace the hadron h in Eq. (6) by a quark state and expand both sides perturbatively, with n, m = 0, 1, 2 indicating the power in α s . While partonic f (n) qv with n = 0, 1, 2 in the MS factorization scheme are known [61], we have to calculate partonic version of F ν(n) qv/q in the MS scheme perturbatively to n = 0, 1, 2 to derive the NNLO matching coefficient K ν(n) .
Some representative Feynman diagrams for F ν(n) qv are shown in Fig. 1 (a , b , c ). The diagram 1(a ) gives the tree level result To calculate F ν qv/q (ξ 2 µ 2 , ω) at high orders, we again use transformation as Eq. (12) to remove the exponential by going to momentum space, and then reduce the loop integrals to MIs by using IBPs. For example, at NLO we have two MIs: and at NNLO we have 21 MIs. The MIs generated from n-loop diagrams for F ν(n) qv/q are functions satisfied where y ≡ p z /q z , d n ≡ −2n − 1, and d n + d nj are the dimensions of MI I (n) j . These MIs can be solved by the differential equations [62] ∂ y J (n) with J (n) j (0; d) serving as boundary conditions. By applying IBPs again, the integrals in boundary conditions can be decomposed into vacuum MIs at n-loop order, which have been calculated in the renormalization procedure. Therefore, J (n) j can be expanded as a Taylor series of y based on the differential equations in Eq. (20).
After carrying out MIs, we can Fourier transform back to position space and the y dependence is changed to dependence on ω. By adding contributions from all diagrams and then multiplying it by UV renormalization factor Z −1 MS , we obtain perturbative results of F ν(n) qv/q (ω, ξ 2 , µ 2 ) with n = 1, 2. We then obtain MS matching coefficients K ν(n) (xω, ξ 2 , µ 2 ) using Eq. (16). As expected, all divergences are canceled and final results of K ν(n) are finite. It verifies the proof of the factorization theorem [13] up to two-loop order. Our one-loop results K ν(1) agree with previous calculations [30][31][32]; terms proportional to n f in two-loop results have been calculated in Ref. [55] using quark mass regulator; while other two-loop results are new. All results up to two-loop order are given in the Appendix.
Using Eq. (6), one can obtain NNLO matching coefficients in other RS by calculating corresponding R RS .

NUMERICAL RESULTS
Based perturbative results calculated above, we present numerical predictions for the valence coordinatespace QCFs by using Eq. (6) with CT18NNLO PDFs as input [63]. We set µ = 2c/|ξ| to minimize logarithms encountered in perturbative calculation when c = 1, and vary c from 1/2 to 2 to estimate theoretical uncertainties due to ambiguity of scale choice. In Fig. 2, we show i 4ω ξ · F vac qv/h (ω, ξ 2 ) as a function of ω with fixed 1/|ξ| = 2 GeV or as a function of 1/|ξ| with fixed ω = 10. In either case it is evident that our numerical results have an improved uncertainty when higher order matching coefficients are used. Especially for small 1/|ξ|, which is the dominant region of lattice data, the NNLO results can reduce theoretical uncertainty by more than a factor of 3 comparing with NLO results.

SUMMARY
We showed that properly renormalized quark correlation functions in position space are good LQCD observables, if ξ 2 Λ 2 QCD is sufficiently small, which are calculable in LQCD and factorizable to PDFs. We discussed the ambiguity and scheme-dependence in defining the multiplicative renormalization constant Z RS , and demonstrated that Z RS defined with the vacuum state is advantageous for carrying out the perturbative calculations of the matching coefficients, especially, at high order in α s . For the first time, we derived a complete NNLO valencequark coefficient functions for QCFs, and demonstrated that the extraction of PDFs from LQCD calculations in terms of QCD factorization approach are in fact at the same rigor as the program to extract PDFs from experimental data.
Our definition of QCFs and method of calculations can be easily generalized to evaluate gluon correlation functions and sea-quark coefficient functions, which provide a rigorous program to extract PDFs from lattice QCD calculations. With multiple "good" LQCD observables, including the current-current correlators (better UV behavior) [6], the quark correlation functions defined here, and a generalization to gluon correlation function in po-sition space, the QCD factorization of these observables provide a complementary revenue for extracting PDFs or other partonic structures of hadrons, as well as a tremendous potential to extract partonic structure of hadrons that could be difficult to do scattering experiments with. Note added: When this work was being finalized, some related preprints appeared [64][65][66]. In Ref. [64] the authors obtained NNLO results for Z MS and R vac , which exactly agree with our results. In Refs. [65,66] the authors obtained matching coefficients for flavor nondiagonal quark to quark channel that starts from twoloop order.
After the manuscript was posted on arXiv, a preprint appeared soon [67], where matching coefficients of valence quark are also calculated to NNLO.