Master integrals for two-loop QCD corrections to quark quasi PDFs

We compute the master integrals for the two-loop QCD corrections to quark quasi parton distribution functions (PDFs) in the large momentum effective theory (LaMET). With a proper canonical basis, we derive the analytical results for the three families of master integrals using the method of differential equations. The final expressions for the master integrals are given in terms of Goncharov polylogarithms. These results allow us to extract the two-loop short-distant matching coefficients between quark quasi and lightcone PDFs in LaMET, and are valuable to improve the determination of the nucleon PDFs from first principles in future.


Introduction
In high energy physics, theoretical predictions for physical observables like cross sections are usually made under QCD factorization scheme, in which the scattering amplitude is split into a perturbative hard-kernel and low-energy matrix elements. The perturbative coefficient arises from the short-distant degrees of freedom, while the long-distant inputs, parton distribution functions (PDFs) and others, describe the longitudinal momentum distribution of unpolarized/polarized partons inside a hadron. The involved partons in the hadron move nearly at the speed of light. Thereby it is extremely difficult to directly calculate the PDFs from first-principles of QCD, i.e. Lattice QCD. Previous attempts in Lattice QCD make use of the operator product expansion and calculate the moments of PDFs. This approach was successful only for the lowest few moments of the lightcone PDFs [1], but studies of higher moments suffer from significantly large noises in the simulation. Thus a complete description of PDFs from first-principles is not established yet.
Recently a breakthrough to circumvent the above problem was proposed in refs. [2,3], and is now formulated as the large momentum effective theory (LaMET). In this framework, instead of calculating the moments of lightcone PDFs one can explore the equal-time correlators, also called quasi observables, on the Lattice. The equal-time correlators have the same infrared structure with lightcone quantities, and more particularly under the large Lorentz boost, the quasi correlators approach the pertinent lightcone observables. Though JHEP10(2020)079 they still have different ultraviolet behaviors, such differences are compensated by the perturbatively calculable coefficients. In the same spirit, other proposals like the "good lattice cross-section" [4,5] and Ioffe-time "pseudo-distributions" [6] are also given in recent years.
Under the LaMET framework, tremendous progress has been made in recent years , and a lot of fairly good results are obtained, many of which are consistent with the phenomenological results extracted from the experiments [65] (see refs. [66,67] for recent reviews). It is anticipated that with the increase of computing resources and the development of various techniques in future, Lattice simulations of quasi distributions will become more and more accurate. On the other hand, most of the available extractions of lightcone PDFs resort to the perturbative kernels at one-loop accuracy. Thus to further reduce theoretical uncertainties, the next-to-next-to-leading order (NNLO) predictions of the matching coefficients are inevitably required. In ref. [63], we have for the first time presented a NNLO calculation for the flavor non-diagonal contributions to the quark quasi distribution functions and explicitly validated the factorization hypothesis at NNLO accuracy. Many calculation details are not included in ref. [63].
Moreover the study of the mathematical properties of Feynman integrals has attracted increasing attention in the last decades, and significant progress was achieved in understanding the analytical behavior of multi-loop Feynman integrals. One of the most powerful tools to evaluate the master integrals analytically is the method of differential equations [68][69][70][71][72]. Along with the recent developments [73][74][75][76], this method has been widely applied to various processes. It has been pointed out in ref. [73] that in a generic multiloop calculation, a suitable basis (canonical basis) of master integrals can be chosen. In this canonical basis, the corresponding differential equations are greatly simplified, and their iterative solutions become straightforward in terms of the dimensional regularization parameter = 4−D 2 . In addition, the determination of the boundary conditions is also greatly simplified.
In this work, we will present the detailed results for the two-loop Feynman integrals of QCD corrections to quark quasi PDFs. With the integration-by-parts (IBP) technique, all the two-loop QCD corrections to quark quasi PDFs are reduced into a set of integrals, called master integrals, which are then calculated using the method of differential equations. In section 2, we will set up the notations and conventions. Section 3 is devoted to the canonical basis of the master integrals, and section 4 details the boundary conditions for the integrals. Section 5 contains the analytical results for the master integrals and some validations through FIESTA are also given. A brief summary is given in the last section.

Notations and conventions
The quasi PDF for the quark q i in a hadron H is defined as: 1) with N being the normalization factor and y ∈ [−∞, ∞] beingẑ-component momentum fractions of the hadron carried by the parton. The q i is the quark field and the W (z, 0) is JHEP10(2020)079 the Wilson line from 0 to z to maintain the gauge invariance. The Γ is the Dirac matrix for the quasi PDF and two popular choices are Γ = γ 0 , γ z for unpolarized PDFs. The gluon quasi PDFs can also be defined similarly.
Studying the quasi PDFs in perturbation theory requires the replacement of hadronic state by partonic state. Then one can calculate the perturbative contributions and extract the short-distant matching coefficients order by order. In the momentum space, the twoloop QCD corrections to quasi PDFs contain a momentum conservation in z direction arising from the Fourier transformation as shown in eq. (2.1). This results in a delta function. In order to use the standard IBP technique, we employ the identity All the involved integrals in the two-loop contributions can be expressed by the three families of integrals, which are generically parameterized as with P 1 = k 2 1 , P 2 = k 2 2 , P 3 = (k 2 − p 1 ) 2 , P 4 = (k 1 + k 2 ) 2 , P 5 = (k 1 + k 2 − p 1 ) 2 , P 6 = n 1 · k 1 + yp z , (2.4) (2.6) The integration measure is defined as and D = 4 − 2 . The momentum in above equations can be parameterized as When p 2 1 ≤ 0 (space-like and lightlike), all the integrals defined above are real.
After some algebraic manipulations and simplifications of the Feynman amplitudes, we encounter a number of tensor integrals. We will use FIRE packages [77][78][79] to perform the IBP reduction, after which all the integrals are reduced to the so-called master integrals. The first family of master integrals contains 36 linearly independent integrals, while the second and third family contain 32 and 28 integrals, respectively. To obtain the analytic results of master integrals, we will apply the method of differential equations and choose the canonical basis.
On the Lattice, the numerical simulation of quasi PDFs is performed in the discretized space-time, and to properly use the simulated results one needs to renormalize the bare quantities. A widely-adopted renormalization approach in the study of quasi PDFs is the so-called regularization-independent momentum subtraction scheme (RI/MOM) [80]. In this scheme, a counterterm with off-shell quarks p 2 1 = −µ 2 R < 0 is introduced. Thus to match the RI/MOM results from Lattice QCD onto the MS PDFs in the continuum theory, one needs to calculate the two-loop diagrams in both cases with p 2 1 = 0 and p 2 1 < 0.
All the rational matrices (M i , N i , L i ) as well as the electronic form of the canonical basis are presented in the Supplementary material.

Boundary conditions
In order to obtain the analytic results for the basis given in the last section, one needs to determine the boundary conditions first. In this section we present the determination of these boundary conditions.

JHEP10(2020)079
can be obtained by directly performing the integration. For the reader's convenience, we show the results for the g 1 1 , g 1 2 , g 1 3 as follows:  The above results are valid in the entire range of y: (−∞ < y < ∞).
Based on the one-loop results for the matching coefficients for the quark quasi PDFs [14,59], we notice that in the (0 < y < 1) region (the so-called physical region which overlaps with the kinematic region for the lightcone PDFs), the analytic results have logarithmic divergences at p 2 1 = 0(z = 1), but are regular at p 2 1 = 4y(y − 1)p 2 z (z = 2y − 1). In the (y < 0, y > 1) regions the perturbative results are regular at p 2 1 = 0(z = 1), but are singular at p 2 1 = 4y(y − 1)p 2 z (z = 2y − 1). It is also interesting to notice that in the space-like region with p 2 1 = −p 2 z (z = 0), the one-loop integrals are also regular. Though the two-loop integrals are much more complicated, we expect these analyticity properties also hold which will be used to determine the boundary conditions for two-loop integrals. The analytic results obtained with this assumption have been validated in the comparison with the results from FIESTA.
For illustration, we consider g 1 7 and g 1 8 , whose differential equations are Since all integrals do not have any singularity at p 2 1 = −p 2 z (z = 0), we find that the above equations imply g 1 7 = 0 at p 2 1 = −p 2 z (z = 0). In the physical region (0 < y < 1), the integrals (g 1 7 , g 1 8 ) are singular at p 2 1 = 0(z = 1) but regular at z = 2y − 1. The above differential equations give 6g In the y > 1 or y < 0 regions, the above integrals are regular at z = 1 and singular at z = 2y − 1. Thus one can obtain Accordingly the boundary condition for {g 1 7 , g 1 8 } is fixed in the entire range of y.

JHEP10(2020)079
Similar with the above examples, all the remaining unknown boundary conditions can be determined from regular conditions at z = {2y − 1, 2y + 1, y, 1, 0}, respectively. We will briefly present the determination of boundary conditions in the following.

The rest boundary conditions for the integrals in the first family
In the first family, the {g 1 6 , g 1 8 , g 1 10 , g 1 12 , g 1 21 , g 1 22 , g 1 27 , g 1 30 } do not vanish at z = 0. The boundary conditions for {g 1 8 , g 1 30 } are determined from the regularity of the corresponding basis at z = 2y − 1 in the (0 < y < 1) region and from the regularity at z = 1 in the (y < 0, y > 1) regions. The boundary conditions for {g 1 12 , g 1 27 } can be similarly obtained. But since these integrals have contributions from the same topology with anti-quark contributions, their boundary conditions are determined from the regularity at z = 2y + 1 in the (−1 < y < 0) region and from the regularity at z = 1 for (y < −1, y > 0).
To derive the boundary conditions as well as analytic results for the g 1 6 and g 1 10 , we consider the region (0 < y < 1) first. In this region the boundary condition for the g 1 6 is determined from the regularity at z = 2y − 1, and then one can obtain the corresponding analytic result straightforwardly. In the y > 0 region, g 1 10 is both singular at z = 2y + 1 and z = 1, and the boundary condition can be determined from the differential equations of g 1 15 at z = y: where the ellipses stand for less singular terms at z = y. Since the p 2 1 = (y 2 − 1)p 2 z (z = y) is space-like, and g 1 15 is regular at p 2 1 = (y 2 − 1)p 2 z , one has Then we can obtain the analytic results for the g 1 10 in the 0 < y < 1 region and the results can be extrapolated to the y > 0 region. In the y < −1 region, the results of g 1 10 can be obtained with the replacement (y → −y − 1) from the ones in the y > 0 region, with an overall minus sign. For the (−1 < y < 0), the g 1 10 is singular at z = 1 but regular at z = 2y + 1, and thus we can determine the boundary condition at z = 2y + 1. This completes the boundary conditions for the g 1 10 in the entire range of y. Results for the g 1 6 in the (y < 0, y > 1) regions can be derived from the ones for g 1 10 in the region (y < −1, y > 0) through the replacement y → y − 1.

JHEP10(2020)079
The base g 1 22 depends on y only, in the region (y < −1, y > 0) it can be determined from the regularity of g 1 26 at z = 1, whose differential equation is given as and the results of g 1 22 in y > 0 region can be extrapolated to the y > −1 region. The results for the g 1 21 are related to the ones for the g 1 22 through y → y − 1.

The rest boundary conditions for the integrals in the second family
For the second family, the boundary conditions for {g 2 9 , g 2 11 , g 2 13 , g 2 15 , g 2 18 , g 2 22 } are to be determined. The boundary conditions for the g 2 11 and g 2 15 as well as full analytic results can be obtained from the ones in the first family through the symmetry relations in eq. (3.11). Results for the g 2 9 can be derived from g 1 22 by reversing the sign for the variable y. The boundary for the {g 2 13 , g 2 18 , g 2 22 } can be determined from the regularity at z = 2y−1 for (0 < y < 1) region and from the regularity at z = 1 for (y < 0, y > 1).

The rest boundary conditions for the integrals in the third family
For the third family, the boundary conditions for {g 3 16 , g 3 18 , g 3 20 , g 3 24 , g 3 27 } need to be determined. Their boundary conditions as well as full analytic results can be obtained from the first and second families using the symmetry relations in eq. (3.11).
The analytic results for all the canonical basis in the entire range of y can be found in the Supplementary material of this work. The explicit expressions of boundary conditions can be obtained directly from the analytic results.

Analytic results and validations
In the Lattice QCD study of quasi PDFs, the RI/MOM scheme is often adopted to renormalize the bare quantities [80]. In this scheme, a counterterm with off-shell quarks is introduced and thus in this work, we will present the analytic results with both off-shell and on-shell quarks.

Results for off-shell quarks
After determining all the boundary conditions, we can obtain the analytic results for all integrals in the entire range of y. We calculate all the integrals up to weight-3, which are required for the involved two-loop corrections. For the first family, we can divide the y range into 4 regions (y < −1, −1 < y < 0, 0 < y < 1, y > 1), while for the second and third family, the results are divided into 3 regions (y < 0, 0 < y < 1, y > 1).

JHEP10(2020)079
For the convenience of the reader, we show the analytic results of some typical integrals at (0 < y < 1) as follows:

JHEP10(2020)079
The results at (p 2 1 = 0, y = 1.35, p z = 1) are Analytic: One can see from the above comparison that the numerical results obtained with the FIESTA perfectly agree with the analytic results.

Discussions and conclusions
In summary, we have presented the detailed results for the two-loop master integrals of NNLO corrections to quark quasi PDFs. Three families of master integrals are derived. Making use of the method of differential equations along with the choice of canonical basis, we obtain the analytic results for all the master integrals in terms of Goncharov polylogarithms and polylogarithms. Our analytic results are in agreement with the numerical results by FIESTA [89,90] package in the entire range of y. These two-loop master integrals are helpful to extract the two-loop matching coefficients between quasi and lightcone PDFs, and accordingly together with the Lattice QCD simulations deepen understanding the lightcone structures inside a hadron. These results have been sued to calculate the NNLO corrections to quasi PDFs distributions for valence quarks [91]. Based on the previous calculation of perturbative kernel at one-loop accuracy, we find that the PDFs for quarks, gluons, meson distribution amplitudes and generalized parton distributions have used some common integrals. Thus from this viewpoint our results in this work with some generalizations are also likely useful in the studies of these quantities. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.