Quasi parton distribution functions at NNLO: flavor non-diagonal quark contributions

We present a next-to-next-to-leading order (NNLO) calculation of the quasi parton distribution functions (Quasi-PDFs) in the large momentum effective theory (LaMET). We focus on the flavor non-diagonal quark-quark channel and demonstrate the LaMET factorization at the NNLO accuracy in the modified minimal subtraction scheme. The matching coefficient between the quasi-PDF and the light-cone PDF is derived. This provides a first step towards a complete NNLO analysis of quasi-PDFs and to better understand the nucleon structures from the first principle of QCD.


INTRODUCTION
Understanding the underlying structure of nucleons from degrees of quarks and gluons has been a long-standing goal in hadron physics. Since deep-inelastic scattering experiments at Stanford Linear Accelerator Center in late 1960's, the proton structure has been explored in various hard scattering processes [1]. The key results involve the parton distribution functions (PDFs), defined as momentum distributions of quarks and gluons in an infinite-momentum hadron. These distribution functions are normally referred as the light-cone PDFs or the collinear PDFs. In high energy experiments at the lepton-hadron and hadron-hadron colliders, the PDFs are also the important ingredients to characterize the structure of hadrons and make predictions for various processes to test the standard model and probe the new physics beyond. Though the scale evolution of PDFs beyond leading order (LO) into next-to-next-to-next-to leading order (NNNLO) have been performed in literatures [2][3][4][5][6], calculating the PDFs and more generally light-cone observables from first principle of quantum chromodynamics (QCD), has been extremely difficult. In the formulation of non-perturbative QCD on a Euclidean lattice, one cannot directly explore time-dependent correlations. Instead, only moments of parton distribution functions, matrix elements of local operators, can be calculated. However, the difficulty in Lattice QCD study grows significantly for higher moments due to technical reasons and thus only limited moments can be extracted to date [7][8][9][10].
An effective theory, called large momentum effective theory (LaMET) [11,12], has been developed to compute various parton distribution functions on Lattice. In this framework, an appropriate static-operator matrix element (quasi-observable) that approaches the parton observable in the infinite momentum limit of the external hadron is constructed. The quasi-observable constructed in this way is usually hadron-momentum-dependent but timeindependent, and thus can be readily computed on the lattice. After the renormalization, the quasi-observable can be used to extract the parton observable through a factorization formula accurate up to power corrections that are suppressed by the hadron momentum. The relevant parton distribution functions calculated in the LaMET are referred as Quasi-PDFs. Great progress has been made in the last few years on both the theoretical understanding of the formalism and the lattice simulations for parton distributions of baryons and mesons, see, for example, some recent reviews in Ref. [13,14].
The factorization arguments of LaMET allow us to carry out order by order perturbative calculations on the matching between the Quasi-PDFs and the light-cone PDFs. This matching is one of the crucial elements in applying LaMET to parton physics. It provides a solid foundation to compute the light-cone PDFs in a systematically controlled way. In some sense, the improvement on the precision of the PDF calculations can only be achieved by combining the advanced lattice simulations for the Quasi-PDFs (toward small lattice spacing, large volume and physical pion mass) and higher order perturbative matching calculations.
Higher order perturbative calculations are also important to demonstrate the factorization in the LaMET explicitly.
In particular, some specific features of the factorization can only be manifest in the non-trivial two-loop calculations.
Quasi-PDFs at one-loop order and the associated matching coefficients has been a subjective of active research since LaMET was proposed in 2013. This includes quark distribution [15][16][17][18], gluon distribution [19][20][21] and many others (See the review [14]). The goal of this paper is to go beyond the one-loop order and perform, for the first time, a two-loop computation of the Quasi-PDF in the LaMET, taking the non-diagonal quark-quark channel as an example. This channel starts at two-loop order, which allows us to demonstrate the factorization in an intuitive method. We also notice that recently, the renormalization of Quasi-PDF operators have been studied at two loop order [22,23]. Together with this result, our paper will provide an important step toward a complete two-loop calculation of Quasi-PDF and the associated matching coefficients. The rest of this paper is organized as follows. In Sec. II, we present our main result of non-diagonal quark-quark splitting in LaMET at two-loop order. We will provide a detailed calculations and demonstrate the factorization in detail. Based on these results, we show the matching coefficients at this order in Sec. III. Since the non-diagonal quark-quark splitting only starts at two-loop order, this presents the leading contribution for this channel. Some numeric results will also be presented to illustrate the behavior of the matching coefficients. We will summarize our paper in Sec. IV.

FACTORIZATION AT TWO-LOOP ORDER
We start with the definitions of the light-cone PDF and Quasi-PDF. For the unpolarized quark light-cone PDF, we have where x = k + /p + is the quark longitudinal momentum fraction and p µ = (p 0 , 0, 0, p z ) is the hadron momentum. Similarly, the quasi-PDF for the unpolarized quark is defined as where z is a spatial direction and we will adopt Γ = γ t with the normalization factor N = p z /p t and use the / p projector.
According to the factorization in the LaMET, we can write down the Quasi-PDFsf q/H (x, p z ) in terms of the light-cone PDFs f q /H (y, µ):f where q , q being the partons in the hadron. Thef q/H (x, p z ) is an equal-time correlation while f q /H (y, µ) is lightcone PDF. Thoughf q/H (x, p z ) and f q /H (y, µ) share the same infrared structure, their ultraviolet behaviors are different, and embedded in the short-distance coefficient C qq . Since the short-distance coefficient is insensitive to the incoming hadrons, in the calculation of C qq one can replace the hadron by the partonic state. In this work we will consider the flavor non-diagonal quark contributions and the hadron state |H is replaced by a quark state |q and we have the condition q = q. We plot the Feynman diagrams for flavor non-diagonal quark distributions in Fig. 1. In our computations below, we will apply the modified minimum subtraction scheme (MS) and dimensional regulation with D = 4 − 2 . Under this scheme, we can write the formula for the flavor non-diagonal quark distribution as where both sides are computed with dimensional regulations and (1/ IR ) n represent the Infrared divergences. At NNLO, the matching scheme is given as Here, we have applied the perturbative expansions with T i being each off q/q , C qq , f q /q . Because of the particular feature of non-diagonal quark-quark splitting, each term at the right hand side of the above equation represents only one contribution. In the first term, q has to be q , so that it only has C (2) q/q . For the second term, q has to be a gluon, and the combination is quark-to-gluon splitting f (1) g/q and gluon-to-quark C (1) qg matching. Finally, q in the third term has to be q, representing non-diagonal quark-quark collinear splitting f (2) q/q . We also know that both f (0) q /q and C (0) q/q are Delta functions. Therefore, the above equation can be simplified as Here, C g/q and f (2) q/q are known in the literature, which are listed in the Appendix for reference. The objective of our calculations is to computef (2) q/q and extract C (2) qq . In the perturbative calculations at this order,f (2) q/q contains only IR divergences, which can be expressed as 1/ IR in the dimensional regulation. According to the factorization theorem, the IR divergences inf (2) q/q will be cancelled by that from the right hand side of Eq. (6). In particular, the 1/ 2 IR term will be cancelled by the last term and the 1/ IR by the second and last term. After these cancellations, we are left with a finite term, which will be the matching coefficient at this order. In the calculations of the Feynman diagrams in Fig. 1, we adopt the covariant gauge and find that the final results are independent of the gauge parameter ξ. We have also performed an independent calculation in the axial gauge, where only subdiagram-(c) and subdiagram-(e) in Fig. 1 contribute. The results are consistent with each other for both covariant and axial gauges. We use FIRE [24] packages to reduce all the involved tensor integrals into a set of integrals called master integrals. There are two independent family of master integrals. Method of differential equations [25,26] are applied to calculate those master integrals. By choosing a set of basis called canonical basis [27], calculation of the differential equations are greatly simplified. We obtain the analytical results for all the canonical basis in all region of x. The analytic results have been cross checked with numerical packages FIESTA [28].
As mentioned above, there is no UV divergence in flavor non-diagonal quark quasidistributions and thus it is not necessary to perform the renormalization in modified minimal subtraction scheme. All soft divergences are also cancelled. The collinear divergences in 0 < x < 1 region contain 1/ 2 IR and 1/ IR : where Γ 1 and Γ 2 are defined as These divergences will be cancelled by two parts in Eq. (6): one from the divergences in the convolution of C qg and f (1) g/q at one-loop order, and f (2) qq at two loop. For the collinear divergences in −1 < x < 0 region, one can obtain from Eq. (7) and do the replacement x → −x and add a prefactor -1. Note that one should do the log(p(x)) → log(p(x) 2 )/2 replacement at first to ovoid to produce the imaginary part. The IR cancellation is similar.
The collinear divergences in nonphysical regions are given as: where These divergences are cancelled by the convolution of C g/q as indicated in Eq. (6). Again, we list the result in the Appendix. One can also see the two regions of x > 1 and x < −1 are related by the symmetry of x → −x and an opposite sign.
The nontrivial cancellation of the IR divergences discussed above is an important demonstration of the LaMET factorization. This also provides a cross check of our final result on the matching coefficient, which will be presented in the next section.
We would like to emphasize a number of points before we close this sections. First, the complete cancellation of the collinear divergence depends on the factorization formula for this channel, see, Eq. (6), including the different terms contributing from the right hand side. Second, it also depends on the exact results of lower order perturbative contributions. For example, the scale dependent term in the one-loop matching C (1) gq (see Appendix) plays a crucial role to demonstrate the complete cancellation in the above equations. This emphasizes the importance of a consistent subtraction scheme in the perturbative calculations of Quasi-PDFs and the matching coefficients. Finally, our example of the non-diagonal quark-quark channel shall provide important guideline for future developments on computing the Quasi-PDFs at two-loop order. We will carry out these calculations in a separate publication.

MATCHING COEFFICIENT AT TWO-LOOP ORDER
The matching coefficient C (2) qq is obtained by expanding both sides of Eq. (6) to O( 0 ) order. Because of all the divergences between them have been cancelled explicitly as shown in the previous section, it is straightforward to carry out the calculations for the finite parts.
First, let us show the result for the region of x > 1. The finite part of f (2) q/q is given bỹ where Γ 1 (x) has been defined in Eq. (11) and h 1 (x) is defined as (13) Expanding all the terms in the right part of Eq. (6) into the order of O( (0) ), the NNLO matching coefficient C (2) qq can be extracted as: where g 1 (x) is defined as Similarly, the finite part of f (2) q/q in 0 < x < 1 is given bỹ where Γ 1 and Γ 2 are defined in Eqs. (8,9) in the last section, respectively, and h 2 (x) is defined as Note that the imaginary parts of terms by log(−x) log 2 (1 + x)/2 and Li 3 (1 + x) from the x → −x transform are cancelled. Then we can obtain matching coefficients in the 0 < x < 1 region where P (1) q/q (x) being the two loop splitting function can be found in Appendix, and g 2 (x) is defined as Matching Coeff .  there is a toy model of f (y) = 6y(1 − y) and adopt the scale µ = p z /y. Therein the left one is for 0.1 < x < 1 and the right one is for 1 < x < 3.
In the end, the NNLO matching coefficient C (2) qq in x < 0 can be obtained by replacing x → −x and adding an overall minus sign. It is interesting to investigate the asymptotic behaviour of the matching coefficients at infinity points. Up to O( ), we have, for example, These will lead to a logarithmic divergence when performing the integration of C g/q (y, IR ) at infinity points. However, they do cancel between the integration at positive and negative infinity points and thus we do not need to add prescriptions to the divergence at the integration of infinity points. For the NNLO matching coefficients, we have From these, one can see C q /q (y, µ) at NNNLO matching. We plot the distributions of matching coefficients of C (2) qq and C (1) qg as a function of momentum fraction x in Fig. 2. Therein the renormalization scale is adopted as µ = p z and the lattice realization of p z is in several GeV currently. From it, C qq has a different shape compared with others. Assuming the parameterization form of light cone PDFs as the simplest one ay b (1 − y) c , we can test the convolution between matching coefficients and light cone PDFs.
So we also plot the convolution of C  qq has double logarithms as Γ 2 log 2 ( µ 2 p z 2 ) in physical region, while single logarithms as Γ 1 log( µ 2 p z 2 ) in non-physical region.

CONCLUSION
In summary, we have presented a next-to-next-to-leading order calculation of the quasi parton distribution functions for the flavor non-diagonal quark contributionsf (2) q/q (x, p z µ ) in −∞ < x < ∞. We have demonstrated the LaMET factorization at this order. The matching coefficient is derived under the modified minimal subtraction scheme.
These results shall be directly employed to investigate the sea quark contributions in both nonsinglet and singlet quark distributions at NNLO. This will stimulate further developments toward a complete calculation of Quasi-PDFs at two-loop order and the associated matching coefficients.
j/k (y, IR ). (23) As mentioned in Sec. II, for the non-diagonal quark-quark splitting at two-loop order, we need to consider the contributions from C (1) qg and f (1) g/q for the second term, C (0) qq and f (2) q /q for the third term. C (0) qq is a Delta function. In the following, we list other terms for the reference. To extract the matching coefficient C qq , we will keep some of the terms up to O( ).
The results of x < 0 can be obtained with the replacement x → −x and an overall minus sign.
To get f (2) q/q , one needs to apply the convolution, to obtain the 1/ 2 IR -term, There is also a 1/ IR -term in f (2) q/q [2,5], The above two contribute to the divergences in the third term of Eq. (6).