Gluon Pseudo-Distributions at Short Distances: Forward Case

We present the results that are necessary in the ongoing lattice calculations of the gluon parton distribution functions (PDFs) within the pseudo-PDF approach. We give a classification of possible two-gluon correlator functions and identify those that contain the invariant amplitude determining the gluon PDF in the light-cone $z^2 \to 0$ limit. One-loop calculations have been performed in the coordinate representation and in an explicitly gauge-invariant form. We made an effort to separate ultraviolet (UV) and infrared (IR) sources of the $\ln (-z^2)$-dependence at short distances $z^2$. The UV terms cancel in the reduced Ioffe-time distribution (ITD), and we obtain the matching relation between the reduced ITD and the light-cone ITD. Using a kernel form, we get a direct connection between lattice data for the reduced ITD and the normalized gluon PDF. We also show that our results may be used for a rather straightforward calculation of the one-loop matching relations for quasi-PDFs.


Introduction
Lattice calculations of parton distribution functions (PDFs) are now a subject of considerable interest and efforts (see Ref. [1] for a recent review and references to extensive literature). Modern efforts aim at the extractions of PDFs f (x) themselves rather than their x N moments. On the lattice, this may be achieved by switching from local operators to current-current correlators [2]. A further idea is to start with equal-time correlators [3] .
X. Ji, in the paper [4] that strongly stimulated further development, made a ground-breaking proposal to consider equal-time versions of nonlocal operators defining PDFs, distribution amplitudes, generalized parton distributions, and transverse momentum dependent distributions. In the case of usual PDFs, the basic concept of Ji's approach is a "parton quasi-distribution" (quasi-PDF) Q(y, p 3 ) [4,5], and PDFs are obtained from the large-momentum p 3 → ∞ limit of quasi-PDFs.
Within the pseudo-PDF approach, the matching relations were derived for non-singlet PDFs [22,23,24,25,15]. The strategy of the lattice extraction of non-singlet GPDs and the pion DA using the pseudo-PDF methods was outlined in a recent paper Ref. [26], where the matching conditions for these cases have been also derived. In the present paper, our main goal is to describe the basic points of the pseudo-PDF approach to extraction of unpolarized gluon PDFs, and also to find one-loop matching conditions.
In the gluon case, the calculation is complicated by strict requirements of gauge invariance. In this situation, a very effective method is provided by the coordinate-representation approach of Ref. [27]. It is based on the background-field method and the heatkernel expansion. It allows, starting with the original gauge-invariant bilocal operator, to find its modification by one-loop corrections. The results are obtained in an explicitly gauge-invariant form.
In this approach, there is no need to specify the nature of matrix element characteristic of a particular parton distribution. This means that one and the same Feynman diagram calculation may be used both for finding matching conditions for PDFs (given by forward matrix elements), and for DA's and GPDs corresponding to non-forward ones (see Ref. [26] for an illustration of how this works for quark operators).
The paper is organized as follows. In Section 2, we analyze the kinematic structure of the matrix elements of the gluonic bilocal operators, and identify those that contain information about the twist-2 gluon PDF.
Next, we discuss one-loop corrections. In Section 3, we analyze the gauge-link self-energy contribution and specific properties of its ultraviolet and short-distance behavior. Our results for the vertex corrections to the gluon link are given in Section 4 in the form that is valid both in forward and non-forward cases. The "box" dia-gram is discussed in Section 5. Since our results in this case are rather lengthy, we present just some of them, and in the forward case only. The gluon self-energy corrections are discussed in Section 6.
The subject of Section 7 is the structure of perturbative evolution of the gluon operators and matching conditions. Section 8 contains a summary of the paper.

Matrix elements
The nucleon spin-averaged matrix elements for operators composed of two-gluon-fields (with all four indices non-contracted) are specified by where [z, 0] is the standard straight-line gauge link in the gluon (adjoint) representation (2. 2) The tensor structures for a decomposition over invariant amplitudes may be built from two available 4-vectors p α , z α and the metric tensor g αβ . Incorporating the antisymmetry of G ρσ with respect to its indices, we have where the amplitudes M are functions of the invariant interval z 2 and the Ioffe time [28] (pz) ≡ −ν (the minus sign is introduced for further convenience).
Since the matrix element should be symmetric with respect to interchange of the fields (which amounts to {µα} ↔ {λβ} and z → −z), the functions M pp , M zz , M gg , M ppzz and M pz − M zp are even functions of ν, The usual light-cone gluon distribution is obtained from g αβ M +α;β+ (z, p), with z taken in the light-cone "minus" direction, z = z − . We have i.e., the PDF is determined by the M pp structure, Thus, we should choose the operators with the sets {µα; λβ} that contain M pp in their parametrization.
Note that it is the density of the momentum G(x) ≡ x f g (x) carried by the gluons rather than their number density f g (x) that is a natural quantity in this definition of the gluon PDF. In the local z − = 0 (or ν = 0) limit, the x-integral gives the fraction of the hadron's plus momentum carried by the gluons. In the absence of gluon-quark transitions, this fraction is conserved, which puts a restriction on the gg-component of the Altarelli-Parisi [29] kernel. Namely, it should have the plus-prescription property when applied to G(x).
The decomposition of these combinations (with summation over i) in the basis of the M structures is (2.10) All of them contain the M pp function defining the gluon distribution, though with different kinematical factors. Unfortunately, none of them is just M pp : they all contain contaminating terms. Moreover, the M 3i;i3 matrix element (proposed originally [4] for extractions of the gluon PDF on the lattice) contains three contaminations, while the others have just one addition. In particular, the matrix element M 0i;i0 has M gg as a contaminating term. It is easy to see that where the summation over both i and j is assumed. Hence, the combination may be used for extraction of the twist-2 function M pp . Combining together matrix elements of different types, one should take into account that, off the light cone, these matrix elements have extra ultraviolet divergences related to presence of the gauge link. Due to the local nature of ultraviolet divergences, each matrix element, for any set of its indices {µα; λβ}, is multiplicatively renormalizable with respect to these divergences [30]. However, choosing different sets of {µα; νβ}, we get, in general, different anomalous dimensions.
Thus, it is not evident a priori which linear combinations of these matrix elements are multiplicatively renormalizable. In Ref. [31], it was established that the combinations represented in Eq. (2.6), namely, M 0i;i0 , M 3i;i3 , M 0i;i3 + M 3i;i0 (and also M 0i;i3 − M 3i;i0 ), with summation over transverse indices i, are each multiplicatively renormalizable at the one-loop level.
Furthermore, the combination G i j G i j (with summation over transverse i, j) equals to 2G 12 G 12 , whose matrix elements are multiplicatively renormalizable. As we will see, it has the same one-loop UV anomalous dimension as M 0i;i0 , hence the combination of Eq. (2.12) is multiplicatively renormalizable at the one-loop level. A possible subject for further studies is to investigate if this is true in higher orders.
The combination g αβ M 3α;3β , containing a covariant summation over α and β, was also found to be multiplicatively renormalizable. It is given by and has the largest number (four) of contaminations. The function g αβ M 0α;0β , also involving a covariant summation, was used in the first attempt [32] of the lattice extraction of the gluon PDF. However, as noted in Ref. [31], it is not multiplicatively renormalizable.
In any theory with a dimensionless coupling constant, the matrix elements M(z, p) contain ∼ ln −z 2 terms corresponding to perturbative (or "DGLAP" for Dokshitzer-Gribov-Lipatov-Altarelli-Parisi [33,29,34]) evolution. One may wonder which combinations have a diagonal DGLAP evolution at one loop.
To answer these questions, we have calculated the modification of the original bilocal operator by oneloop gluon exchanges.

Link self-energy contribution and ultraviolet divergences
The simplest diagram corresponds to the self-energy correction for the gauge link (see Fig. 1). Its calculation is the same as in case of the quark bilocal operators (see, e.g., Ref. [23]). At one loop, one should just the change the color factor C F → C A . Thus, we have where D c µν (z) = g µν /4π 2 z 2 is the Feynman-gauge gluon propagator in the coordinate representation. The resulting integrals over the link parameters t 1 , t 2 diverge when t 1 ∼ t 2 , i.e., when the endpoints t 1 z and t 2 z of the gluon propagator are close to each other. So, one may suspect that this divergence has an ultraviolet origin. To see that this is the case, we use the dimensional regularization (DR) [35] in the UV region, switching to d dimensions. As a result, the gluon propagator in the coordinate space acquires an extra factor (−z 2 ) 2−d/2 . This results in an extra (t 2 − t 1 ) 2−d/2 factor in Eq. (3.2), and the integral there converges for sufficiently small d.
To preserve gauge invariance, our calculations were made using massless gluons and the dimensional regularization. However, in the case of the link self-energy diagram, the use of DR (which is basically just a mathematical trick) is rather misleading in a couple of points.
The relevant subtleties may be illustrated by using the Polyakov prescription 1/z 2 → 1/(z 2 − a 2 ) for the gluon propagator in the coordinate representation [36] (see also Refs. [37,23]). It softens the gluon propagator at intervals −z 2 a 2 , and eliminates its singularity at z 2 = 0. In this respect, it is similar to the UV regularization produced by a finite lattice spacing, and gives The regularized integral vanishes on the light cone z 2 = 0 and converges for spacelike z. Taking z = z 3 and calculating the integrals gives [37,23] The result contains a linear ∼ 1/a divergence that is missed if one uses the DR. Furthermore, for a fixed a and small z 3 it behaves like z 2 3 /a 2 , i.e., Γ Σ (z, a) vanishes for z 3 = 0, as expected: there is no link if z 3 = 0. It also vanishes on the light cone z 2 = 0.
The fact that Γ Σ (z 3 = 0, a) = 0 means that, for a fixed a, this term gives no corrections to the local limit of the G µα (z) [z, 0] G λβ (0) operator, e.g., to the energymomentum tensor (EMT). Since the matrix element of the EMT gives the fraction of the hadron momentum carried by the gluons, the link self-energy correction does not change this fraction. This is a natural phenomenon in the absence of the gluon-quark transitions.
However, if one formally takes the a → 0 limit for a fixed z 3 in Eq. (3.4), then ln 1 + z 2 3 /a 2 converts into 3 the expression ln z 2 3 /a 2 singular for z 3 = 0. Similarly, using the DR, one faces an outcome proportional to where µ UV is the scale accompanying this UV dimensional regularization. Again, the starting expression vanishes for z 2 = 0, but renormalizing it by a subtraction of the 1/ε UV pole, one may apparently conclude that, in addition to the UV divergence, this diagram contains a singularity on the light cone z 2 = 0. For this reason, in our DR results we will explicitly separate the z 2 -dependence induced by the UV singular terms (that actually vanish on the light cone) and that present in the DGLAP-evolution logarithms ln −z 2 µ 2 IR , where µ IR is the scale associated with the DR regularization of the collinear singularities.
The main difference is that if, instead of DR, one regularizes collinear singularities by using a physical IR cut-off Λ (like nonzero gluon virtuality or gluon mass), the one-loop result, proportional to the modified Bessel function K 0 ( √ −z 2 Λ 2 ), remains singular for z 2 = 0, unlike the UV-induced logarithm ln 1 − z 2 /a 2 .
In the case of the link self-energy diagram, we have UV singularities only. Its correction to the G µα (z)G λβ (0) operator is given by

Vertex contributions
There are also vertex diagrams involving gluons that connect the gauge link with the gluon lines, see Fig. 2.
We use the method of calculation described in Ref. [27]. It is based on the background-field technique, with the gluon propagator taken in the "background-Feynman" (bF) gauge [27]. It should be noted that the three-gluon vertex in the bF gauge is different from the usual Yang-Mills vertex (see e.g. [38]). Therefore, the results obtained for separate diagrams in the bF gauge differ from those obtained in the usual Feynman gauge and only the sum of all diagrams must be the same.

UV divergent term
Clearly, the gluon exchange produces a correction just to one of the fields in the G µα (z)G λβ (0) operator, while another remains intact. In particular, the diagram 2a changes G µα (z) into the sum of two terms. One of them contains UV divergences, while the other one is UV finite. The UV-divergent term is given by where G zσ ≡ z ρ G ρσ andū ≡ 1 − u. The overall d-dependent factor here is finite for d = 4, but the u-integral diverges at the lower limit. Thus, just like in the case of the link self-energy diagram, the divergence appears in the integral over a dimensionless parameter t specifying the location of the endpoint of the gluon line on the gauge link. The divergence disappears if one uses the UV regularization by taking d = 4 − 2ε UV , which converts it into a pole at ε UV = 0. Since the UV divergence comes from the u → 0 integration, we can isolate it by takingū = 1 in the gluonic field, which gives The remainder is given by where the plus-prescription at u = 0 is defined as  At first sight, the field G µα (z) = z α G zµ (z) − z µ G zα (z) accompanying the UV pole in Eq. (4.2) does not look like the field G µα (z) in the original operator. Thus, one may worry that we are not dealing here with a multiplicative UV renormalization. So, let us perform an explicit check for our particular case when z = {0, 0, 0, z 3 }.
To begin with, we see that G µα (z) = 0 when both µ and α are transverse indices i, j. This corresponds 4 to a multiplicative renormalization with the anomalous dimension (AD) equal to zero.
Thus, for all the cases, G µα (z) is a multiple of G µα (z). Namely, when one of the indices equals 3, we have a nontrivial anomalous dimension, since G 3α (z) = −G α3 (z) = z 2 3 G 3α (z). In all other cases, we have a trivial (vanishing) AD, since G i j (z) = 0 and G 0i (z) = 0. As mentioned, the link self-energy diagram has both linear and logarithmic UV divergences, while the vertex diagrams have just logarithmic UV divergences. Adding the logarithmic UV divergence coming from the link self-energy to the UV divergences of the vertex diagrams, we find, in particular, that the matrix elements M 0i;i0 and M i j;i j have the logarithmic AD due to the link self-energy diagram only. Call it γ. Comparing overall factors in Eqs. (3.6) and (4.2), we conclude that M 3i;i3 has the logarithmic AD equal to 2γ and matrix elements M 0i;i3 ± M 3i;i0 have the logarithmic AD equal to 3 2 γ. In addition, all of these structures acquire at one loop the same factor due to the linear UV singularity.

Evolution term
Our calculations show that the second, UV finite term from the diagram 2a is given by Note that the gluonic operator in Eq. (4.5) has the same tensor structure as the original operator G µα (z)G βν (0) differing from it just by rescaling z →ūz. There is no mixing with operators of a different type. The uintegral in this case does not diverge for d = 4, but the overall Γ(d/2 − 2) factor has a pole 1/(d − 4).
Formally, there is also a pole 1/(d−3), corresponding to a linear UV divergence. However, the singularity for d = 3 is eliminated by the u 3−d − 1 combination in the integrand. One may say that the linear divergences present in "u 3−d " and "−1" parts cancel each other.
In the calculation of Refs. [31,18] performed using the usual Feynman gauge, the linear singularities cancel between contributions of two different diagrams shown in Fig. 1 of Ref. [18]. In our calculation, based on the bF gauge, the sum of these diagrams is represented by just one vertex diagram, so that the cancellation occurs inside the contribution (4.5) of that diagram.
The remaining 1/(d − 4) pole corresponds to a collinear divergence developed because all the propagators correspond to massless particles. Taking a nonzero gluon mass λ, one would get a finite result containing K 0 ( √ −z 2 λ) (see, e.g., Ref. [23] for a discussion of the quark vertex diagram in a similar context).
is only finite as far as z 2 is finite. The IR cut-off does not eliminate the logarithmic singularity ln −z 2 λ 2 that K 0 ( √ −z 2 λ) has on the light cone. In the z = z 3 case, z 2 3 works like an ultraviolet cut-off for this singularity. This may be contrasted with the UV divergent contributions, where the UV cut-off is provided by the Polyakov regularization parameter a (or lattice spacing a L ) while z 2 3 appears on the IR side of the relevant logarithm ln z 2 3 /a 2 .

Box diagram
There is also a contribution given by the diagram in Fig. 3 containing a gluon exchange between two gluon lines. This diagram does not have UV divergences, but it has DGLAP ln z 2 3 contributions. In contrast to the vertex diagrams, the original G µα (z)G νβ (0) operator generates now a mixture of bilocal operators corresponding to various projections of G µα (ūz)G νβ (0) onto the structures built from vectors p, z and the metric tensor g.
In particular, in the case of the original p| G 0i (z)G 0i (0) |p matrix element, the box diagram contribution is expressed through matrix elements of p| G 0i (uz)G 0i (0) |p , p| G 3i (uz)G 3i (0) |p , p| G 30 (uz)G 30 (0) |p and p| G i j (uz)G i j (0) |p types. All these matrix elements also appear in the box diagram if one starts with the p| G 3i (z)G 3i (0) |p matrix element. Thus, in both cases we have a complicated mixing of different types of operators.
The situation is simpler for matrix elements Namely, for M + 03 (z, p) (or M − 03 (z, p)) combination, the box diagram contribution is expressed through M + 03 (uz, p) (or M − 03 (uz, p)) only. However, does not contain the twist-2 function M pp , and is of no interest. For M + 03 (z, p), the box contribution is given by Here, the Γ(d/2 − 2) terms are singular for d = 4, which results in ln −z 2 terms generating the DGLAP evolution. The Γ(d/2 − 1) terms are singular for d = 2, which corresponds to the fact that the gluon propagator in two dimensions has a logarithmic ln −z 2 behavior in the coordinate space. For d = 4, these terms are finite. Note that, unlike the vertex part, the box contribution does not have the plus-prescription form.

Gluon self-energy diagrams
One may expect that the plus-prescription form would appear after the addition of the gluon self-energy diagrams, one of which is shown in Fig. 4a. These diagrams have both the UV and collinear divergences. On the lattice, the UV divergence is regularized by the lattice spacing. In a continuum theory, one may use the Polyakov prescription 1/z 2 → 1/(z 2 − a 2 ) for the gluon propagator. The collinear divergences may be regularized by taking a finite gluon mass λ. The result is a ln a 2 λ 2 contribution. However, it does not have the z-dependence, and apparently cannot help one to build the plus-prescription form for the ln z 2 3 part of the box contribution.
A possible way out is to represent ln a 2 λ 2 as the difference ln z 2 3 λ 2 − ln z 2 3 /a 2 of the evolution-type logarithm ln z 2 3 λ 2 and a UV-type logarithm ln z 2 3 /a 2 . The latter can be added to the UV divergences of the diagrams 1 and 2 corresponding to link self-energy and vertex corrections. The ln z 2 3 λ 2 part is then added to the evolution kernel.
To be on safe side with gauge invariance, we use the dimensional regularization. Then the analog of the ln a 2 λ 2 logarithm is a pole 1/(2 − d/2) sometimes written as 1/ UV − 1/ IR . For our purposes, it is more convenient to symbolically write it in a form similar to ln a 2 λ 2 . Changing λ → µ IR and a → 1/µ UV we get ln µ 2 IR /µ 2 UV , and then split this into the difference ln z 2 3 µ 2 IR − ln z 2 3 µ 2 UV . We should also take into account the diagrams (one of them is shown in Fig. 4b) with an extra gluon line going out of the link-gluon vertex. The combined contribution of the Fig. 4 diagrams and their left-leg analogs is given by where β 0 = 11N c /3 in gluodynamics, so that the terms in the square bracket combine into 1/6. As discussed above, we will treat 1/(2 − d/2) as ln z 2 3 µ 2 IR − ln z 2 3 µ 2 UV .

When DGLAP is diagonal in pure gluodynamics
The M + 03 ≡ M 0i;i3 + M 3i;i0 combination defined by Eq. (5.1) contains the twist-2 amplitude M pp , though with a higher-twist admixture M zp + M pz . In the local limit, the relevant operator is proportional to the 3rd component of the Poynting vector As already mentioned, the box part of the one-loop correction to the matrix element M + 03 (z 3 , p) in pure gluodynamics has a simple DGLAP structure 1 (5.3). Combining all the gluon one-loop corrections to it, we get, in the MS scheme, The first line here comes from the UV-singular contributions. It contains the δ(ū) factor which reflects the local nature of the UV divergences and converts M + 03 (uz, p) into M + 03 (z, p). The second line contains the Altarelli-Parisi (AP) kernel It has the plus-prescription structure reflecting the fact that, in the local limit, M pp (z, p) is proportional to the matrix element of the gluon energy-momentum tensor. From now on, "+" means the plus-prescription at 1. The third line contains z 3 -independent terms coming from the vertex diagrams (these have the plusprescription form) and from the box diagram. The latter may be written as a sum of the term 2 ūu + 2 3ū 3 + that has the plus-prescription form and the term 2 3 δ(ū) that may be combined with the UV terms.

Reduced Ioffe-time distribution
The combination M pz + M zp is an odd function of ν = z 3 p 3 . Writing it as 2z 3 p 3 m + zp (ν, z 2 3 ), we have Dividing out the kinematical factor 4p 3 p 0 , we deal with which is a function of ν and z 2 3 . Now, just like in the quark case considered in Refs. [10,12], we can introduce the reduced Ioffe-time distribution Since M pp (ν, z 2 3 ) is obtained from the multiplicatively renormalizable combination M + 03 , the UV divergent Z(z 2 3 µ 2 UV ) factors generated by the link-related and gluon self-energy diagrams cancel in the ratio (7.6). As a result, the small-z 2 3 dependence of the reduced pseudo-ITD M(ν, z 2 3 ) comes from the logarithmic DGLAP evolution effects only. Moreover, M pp (0, z 2 3 ) has no DGLAP logarithmic dependence on z 2 3 , because of the plus-prescription nature of the AP kernel B gg (u).
Thus, neglecting O(z 2 3 ) terms, we conclude that, in pure gluodynamics, M(ν, z 2 3 ) satisfies the evolution equation with respect to z 2 3 . This relation is modified when gluon-quark transitions are present.

Gluon-quark mixing
In the MS scheme, the contribution to M + 03 from the gluon-quark diagram shown in Fig. 5 is given by where O q (z 3 ) is a singlet combination of quark fields, with f numerating quark flavors. Note that O q (z 3 ) vanishes for we see that the lowest term is proportional to the quark part of the energy-momentum tensor. This term is accompanied by the z 3 factor which cancels the overall 1/z 3 factor in Eq.

Matching relations
A disadvantage of the M + 03 (z 3 , p) combination is that it vanishes when p 3 = 0 (see Eq. (7.4)). Thus, to extract M pp (ν, z 2 3 ) for ν = 0, one should make measurements of M + 03 (z 3 , p) for a few low values of p 3 , divide p 3 out, and extrapolate the results to p 3 = 0. This procedure leads to additional systematic uncertainties.
Fortunately, the combination M 0i;i0 − M i j;i j = 2p 2 0 M pp of Eq. (2.12), being proportional to p 2 0 , does not have this problem. Furthermore, it gives the twist-2 amplitude M pp without contaminations. The amplitude M pp (ν, z 2 3 ) obtained in this way may be used to form the reduced pseudo-ITD M(ν, z 2 3 ), as in Eq. (7.6). Using the results of our calculations for the oneloop corrections to M 0i;i0 and M i j;i j , and keeping just the M pp term in the correction (while skipping the "higher twist" terms M zz , M zp , M pz , M ppzz ) we obtain the matching relation dw I S (wν, µ 2 ) − I S (0, µ 2 ) B gq (w) (7.17) between the "lattice function" M(ν, z 2 3 ) and the lightcone ITDs I g (ν, µ 2 ) and I S (ν, µ 2 ). The first of them is related to the gluon PDF f g (x, µ 2 ) by I g (ν, µ 2 ) = 1 2 1 −1 dx e ixν x f g (x, µ 2 ) . (7.18) Since x f g (x, µ 2 ) is an even function of x, the real part of I g (ν, µ 2 ) is given by the cosine transform of x f g (x, µ 2 ), while its imaginary part vanishes. The factor I g (0, µ 2 ) has the meaning of the fraction of the hadron momentum carried by the gluons, I g (0, µ 2 ) = x µ 2 .
Thus, Eq. (7.17) allows to extract the shape of the gluon distribution. Its normalization, i.e., the value of x µ 2 should be found by an independent lattice calculation, similar to that performed in Ref. [39]. The singlet quark function I S (wν, µ 2 ) that appears in the O(α s ) correction should be also calculated (or estimated) independently.
The important point is that the R(y, z 2 3 µ 2 ) kernels are given by explicit perturbatively calculable expressions. Using them and Eq. (7.19) one may directly relate M(ν, z 2 3 ) and the light-cone PDFs. Then, assuming some parameterizations for the f g (x, µ 2 ) and f S (x, µ 2 ) distributions, one can fit their parameters and α s from the lattice data for M(ν, z 2 3 ) using Eqs. (7.19), (7.20). This procedure is essentially the same as that used in the "good lattice cross sections" approach [6,7].