Renormalization of parton quasi-distributions beyond the leading order: spacelike vs. timelike

We argue that the renormalization factors for nonlocal quark-antiquark and gluon operators at space-like and time-like separations connected by a Wilson line coincide to all orders in perturbation theory. We calculate the anomalous dimensions and renormalization constants of quark-antiquark and gluon operators to three- and two-loop accuracy, respectively, and also compute vacuum expectation values of these operators to three-loop accuracy.


Introduction
Studies of non-local gauge-invariant operators containing segments of Wilson lines have a long history. They emerged in connection with the attempts to reformulate gauge theories in terms of path-ordered gauge factors (Wilson lines), in particular within the loop-space formalism by Makeenko and Migdal [1]. The study of the renormalization of Wilson lines has been initiated by Polyakov [2], Gerwais and Neveau [3], and continued by several authors [4][5][6][7]; see Ref. [8] for a review. The one-dimensional auxiliary-field formalism introduced in this context in Refs. [3,6] enables the application of the usual language of correlation functions of local operators and is important both conceptually and at a technical level, as a basis for multiloop calculations. Subsequent applications of these methods have been to the study of infrared singularities in Feynman amplitudes, starting from the work of Refs. [9,10], heavy-quark effective theory (HQET) [11][12][13][14][15] and transverse-momentum-dependent (TMD) factorization [16].
Recently, there has been renewed interest in the study of matrix elements of nonlocal off-light-cone operators of the type Q(z) =q(zv) Γ[zv, 0] q(0) , (1.1) where q(x) is a quark field, F µν (x) is the gluon field strength tensor, Γ is a certain Dirac structure, v µ is an auxiliary four-vector and z is a real number. In addition, [zv, 0] is a straight-line-ordered Wilson line connecting the two fields, which we assume to be taken in the proper representation of the gauge group, fundamental for quarks and adjoint for gluons. Matrix elements of such operators acting on hadron states with large momenta are often referred to as parton quasi-distributions (qPDFs) [17] or pseudo-distributions [18]. They can be factorized in terms of parton distribution functions (PDFs) [19] and, at the same time, are accessible in lattice calculations if the quark-antiquark separation is chosen to be space-like, v 2 ≡ v µ v µ < 0. It was suggested that PDFs can be constrained in this way [17], and this possibility is being intensively explored, see, e.g., Refs. [20,21] for reviews. The rationale for using matrix elements of the operators in Eqs. (1.1) and (1.2) in these studies is that they are "cheaper" to compute on the lattice as compared to other Euclidean observables with similar factorization properties; see, e.g., Refs. [22][23][24]. This advantage comes at the price that the renormalization of the non-local operators in Eqs. (1.1) and (1.2) is nontrivial and requires special attention [19,[25][26][27][28][29][30][31][32][33].
In this paper, we address the question whether computational methods familiar from HQET can be applied to the calculation of the renormalization constants (RCs) of the operators in Eq. (1.1), alias for qPDFs, in high orders. The difference is that, in HQET, the heavy-quark velocity v µ is time-like, v 2 > 0, while, for qPDF studies, it is space-like, v 2 < 0. We argue that the change of sign has no effect on the renormalization and confirm this result by an explicit calculation to three-loop accuracy. This result is also relevant in the context of TMD factorization, where Wilson lines are shifted off the light cone to regularize rapidity divergences in TMD operators [16]. Our statement is that the anomalous dimensions (ADs) and RCs do not depend on the direction, space-like or time-like. To avoid confusion, in this work, we imply using dimensional regularization. In renormalization schemes with an explicit regularization scale, the Wilson line in Eq. (1.3) suffers from an additional linear ultraviolet divergence [2], which has to be removed. This can be done by mass renormalization, similarly to the introduction of the residual mass term in HQET, or, alternatively, by considering a suitable ratio of matrix elements involving the same operator [34,35]. Having in mind the second approach, we calculate in this work the vacuum expectation values (VEVs) of the operators in Eqs. (1.1) and (1.2) to three-loop accuracy for space-like and time-like choices of the auxiliary vector v µ . VEVs of nonlocal gluon operators are also of interest for studies of the QCD vacuum structure; see Ref. [36] for a review and further references. The gluon case is also interesting as the generic gluon operator as defined in Eq. (1.2) is not renormalized multiplicatively. The calculation of its VEV allows one to obtain the renormalization constants avoiding the necessity to consider mixing with non-gauge-invariant operators. In this way, we verify the mixing pattern found in Refs. [37,38] and also calculate the two-loop mixing matrix, which is another new result.
This paper is organized as follows. In Section 2, we recall the standard argument that any Green function in QCD with the insertion of the non-local operator in Eq. (1.1) may be obtained from the correlation function of two appropriate local "heavy-light" operators and verify this diagrammatically. In Section 3, we introduce the formalism to be used here, explain the reduction of Feynman diagrams to master integrals and discuss the relation of Green functions between the time-like and space-like regions. In Section 4, we explain how the relationship between time-like and space-like Green functions translates from momentum space to position space. In Section 5, we present our analytic and numerical results for the correlators of interest here, both in momentum and position space. Section 6 contains our conclusions. For the reader's convenience, we list the analytic results through three loops for the ADs and RCs entering our analysis in Appendices A and B, respectively.

Effective field theory
As is well known, the interaction of a particle propagating along a classical path in the background gauge field reduces to the path-ordered phase factor along its trajectory. Such an auxiliary classical particle can be simulated by supplementing the QCD Lagrangian L QCD by an extra term, where h v is a (complex) scalar field in either fundamental or adjoint representation of the gauge group, D = ∂ µ − igA a µ T a is the covariant derivative and T a are the SU(3) generators in the appropriate representation. The Lagrangian in Eq. (2.1) for v 2 > 0 is, essentially, the standard HQET Lagrangian [39], apart from our choice of h v (x) as a scalar. The case v 2 < 0 is of interest in connection with qPDFs [27,31]. Without loss of generality, we can assume v 0 > 0 and the usual causal boundary conditions for the "heavy" field h v . Its free propagator reads /v 2 and i, j are color indices. Adding the interactions with the gluon field, one obtains [8] so that any Green function in QCD with the insertion of the non-local operator in Eq. (1.1) can equivalently be obtained from the correlation function of two local "heavylight" operators. E.g. for quarks, of the path-ordered exponential in Eq. (1.3), one obtains where we have used the shorthand notation A v = v µ A µ . Notice that the factor 1/2! that naively appears in the third term of the expansion of the exponential function is replaced by the ordering of the integration regions. On the other hand, one can easily write the corresponding expression using the Feynman rules of the effective theory in Eq. (2.1), In this way, the 1/2! factor in the last term is compensated by a different mechanism: there are two equivalent ways to couple the four h v fields in the Lagrangian insertions . Notice that disconnected diagrams are not counted. The two expressions in Eqs. (2.5) and (2.6) are obviously equal up to an overall factor, which is nothing but the free propagator of the "heavy" field, S h (zv).

Generalities
We define the renormalization factors for generic composite operators O as where O B is the corresponding bare operator and the bare fields are related to the renormalized ones as The bare QCD coupling constant g B is expressed as where µ is the 't Hooft mass and = (4 − D)/2, with D being the running spacetime dimension within the method of dimensional regularization [40][41][42]. Within the modified minimal-subtraction (MS) scheme, every RC is independent of dimensional parameters (masses and momenta) and can be represented as where a = g 2 /(16π 2 ). Given a RC Z(a), the corresponding AD is defined as Here, the sign depends on the way how the considered bare quantity, O B , is related to its renormalized counterpart, O. With our convention (3.1), the renormalization group (RG) equation for O assumes the form Traditionally, in renormalizing coupling constants, wave functions and masses, the defining relation (3.1) is written with Z O on the left-hand side, e.g., As a result, the ADs γ 2 and γ h , and the β function (see below) are defined with the minus sign version of Eq. (3.5), but the ADs of composite operators (see below) with the plus sign. Customarily, one also defines Z a = Z 2 g and refers to the corresponding AD as the QCD β function, β n a n . (3.7) The RCs Z 2 and Z a serve to renormalize the standard QCD Lagrangian and have been well known through three loops for a long time [43][44][45]. The RC Z h is known to the same accuracy for the time-like case v 2 > 0 [13] and will be used by us as convenient reference point. For the reader's convenience, we quote the RCs as well as the corresponding ADs in Appendices B and A, respectively. Due to an overall δ (D−1) (x ⊥ ) factor in position space, the self-energy and hence also the propagator of the "heavy" field h v can only depend on the projection of its momentum onto v µ , for which we use the notation The heavy-field self-energy is defined, as usual, as the sum of one-particle-irreducible amputated diagrams, 9) and the full propagator is then given bỹ Notice that the full propagator and the self-energy in Landau gauge satisfy the RG equations, The renormalization factors Z Q of the "heavy-light" operators in Eq. (3.13) are calculated from the corresponding vertex functions Γ(p, ω) as explained, e.g., in Ref. [46]; the general formalism is the same for v 2 > 0 and v 2 < 0.

Renormalization constants and VEVs of qPDF operators
Let us first consider quark operators. Let where i is the spinor index, which we do not show in what follows. One can show that this operator is multiplicatively renormalized. The corresponding RC Z Q , was calculated to three-loop accuracy for the time-like case v 2 > 0 [13] and is the same The corresponding ADs are obviously related as The equality (3.15) implies that Z Q does not depend on the Dirac structure Γ in the definition of the operator Q(z), Eq. (1.1). We have checked this relation at the two-loop level by explicit calculation of the VEV of Q(z) to three-loop accuracy; see Section 5. This VEV is non-vanishing in perturbation theory only for the Dirac structure Γ = γ µ , in which case it follows from Lorentz invariance that 0|Q µ |0 ∝ v µ . Thus, it is sufficient to consider 17) or, equivalently, the momentum space correlation function, where the current Q(x) is defined in Eq. (3.13). It is easy to see that where the ± sign corresponds to the choice v 2 = ±1. Notice that Π(z > 0) does not depend on the "contact" terms inΠ(p) of the form const · ω 2 , which suffer from an extra UV divergence coming from the integration region around x = 0 in Eq. (3.18). Thus, in momentum space, it is natural to consider instead ofΠ(p) an analog of the Adler function, namelyD The corresponding RG equation for D(ω) reads The gluon case is more involved because, in the presence of the external four-vector v µ , the components parallel and transverse to v µ can be renormalized differently. Introducing the corresponding projection operators, we can define two multiplicatively renormalizable gauge-invariant local operators in effective theory with adjoint "heavy" scalars as † with the RCs Z ⊥ and Z ⊥⊥ , We will denote the corresponding ADs as γ ⊥ and γ ⊥⊥ respectively. The correlation functions of these operators have the form whereḠ ⊥⊥ αβ = g ⊥ µα g ⊥ νβ gF αβh v , etc., and are renormalized by squares of the corresponding RCs in Eq. (3.24). In this notation, a generic two-gluon vacuum correlation function related to the qPDF operator in Eq. (1.2) takes the form For future reference, we introduce the corresponding "Adler" functions: The functionsD ⊥⊥ (ω) andD ⊥ (ω) satisfy the standard RG equations like Eq. (3.21) with the ADs (2 γ ⊥⊥ ) and (2 γ ⊥ ) respectively. In the qPDF literature, following Ref. [37,38], one usually introduces a different operator basis, so that, obviously, The operators J 1 and J 2 are renormalized by a triangular 2 × 2 mixing matrix Z ik such that Z 21 = 0 and Z 12 − Z 22 + Z 11 = 0 to all orders. These relations imply that Notice that we ignore mixing with gauge-noninvariant operators as they do not contribute to gauge-invariant observables. With the main definitions at hand, we proceed to describe the calculation procedure.

Reduction
The generation of the Feynman diagrams and their reduction to master integrals (MIs) have been done in the standard way, using the programs QGRAF [47] and FIRE6 [48], respectively. ‡ The reduction typically results in a sum of MIs with coefficients being rational functions of the space-time dimension D. In addition, every MI is multiplied by a factor of the form and color factors like C F , C A , etc. We have checked that our results for the time-like case v 2 = 1 are in full agreement with Ref. [13]. In our calculation, we have used the values of the relevant MIs given in Ref. [13].

Master integrals: time-like versus space-like
The MIs for v 2 > 0, which we refer to as time-like, are analytic functions of ω. They are real for ω < 0 and have a branch cut at ω > 0. On dimensional grounds, each MI has the form where M ( ; n) is a real function of the space-time dimension D, d M is the dimension of the MI M , L stands for the number of loops and the sums over l and h count the indices n l and n h of all usual QCD ("light") and special ("heavy") lines of the integral. The argument n stands for the collection of all indices. We assume that every MI is of scalar type, so that the corresponding integrand is given by a product of denominators involving "light" (massless) and "heavy" propagators, possibly raised to certain (integer) powers (indices). The reduction to scalar MIs is certainly possible at the three-loop level; see, e.g., Refs. [13,51]. The restriction to the normalization v 2 = 1 can easily be relaxed. Indeed, the "heavy" propagator in Eq. (2.2) is a homogeneous function w.r.t. the rescaling v µ → λv µ . Thus, we have M (λω, The result for ω > 0 is obtained by analytic continuation. In this way, the MI acquires an imaginary part according to the usual causal prescription ω → ω + i0, so that 4ω 2 → (−2ω − i0) 2 . In order to calculate a MI for the space-like case v 2 < 0, it is useful to start from the so-called α representation of the time-like MI, § where Σn = α n α ≡ l n l + h n h , Σz = α z α , and U and F are homogeneous polynomials of degree L and L + 1 in the integration parameters z α , respectively. The function U does not depend on kinematic invariants, whereas the function F can be written as with polynomials T p and T v that only depend on the parameters z α and are defined to be positive in the integration region of Eq. (3.35). Notice that, in Eq. (3.35), we do not assume that ω < 0, so that, for v 2 > 0, the MI acquires an imaginary part at ω > 0 according to the Feynman prescription F → F − i0. The crucial observation is that, if one simultaneously changes the signs of v 2 and ω, the MI receives an overall phase factor, as, by definition, Thus, we obtain given by the sum of MIs multiplied by extra kinematic factors (v 2 ) j 1 (−2ω) j 2 , where j 1 and j 2 are integers that satisfy the obvious relations It is easy to see that so that, upon this multiplication, the mass and v dimensions of a particular MI are substituted by those of the Green function in question and are the same for the contributions of all MIs and for all Feynman diagrams. As a consequence, going over from v 2 > 0, ω < 0 to v 2 < 0, ω > 0, the Green function acquires an overall phase factor, , which is the final result. Thus, we conclude the following: • A generic Green function at v 2 < 0 can be obtained from the result at v 2 > 0 by the (possible) global sign change (−1) (d G 0 +d G v )/2 , which is the same to all orders of perturbation theory, and the formal substitution where we assume |v 2 | = 1 and the Feynman causal prescription.
• Since, in minimal schemes, the RCs and ADs neither depend on the global sign nor on the value of the external momentum, they are not affected by the analytic continuation in v 2 and are the same for time-like (v 2 = 1) and space-like (v 2 = −1) kinematics.

Transition to position space
Without loss of generality, we may assume v µ = (1, 0, 0, 0) and v µ = (0, 0, 0, 1) for the time-like and space-like cases, respectively. For this choice, the variable ω is the energy, ω = p 0 , for the time-like case and the z component of the momentum up to a minus sign, ω = −p z , for the space-like case. The relation between generic correlation functions established in the previous section, therefore, connects a v 2 = 1 Green function for negative energy with the corresponding v 2 = −1 Green function with negative momentum in z direction. The corresponding position-space variables for these cases are obviously the separation in time t and distance z of the quark and the antiquark in the operator in Eq. (1.1). In the time-like case, the transition is performed with the help of the generic formula where we assume n to be integer. Thus, renormalized time dependent correlation functions are expressed in terms of the following combination where Euler's constant γ E appears naturally due to a universal factor, with (1 + 2L ) (n−1) being the Pochhammer symbol. There is no γ E in momentum space results. The analogue of Eq. (4.1) for the space-like case is Comparing Eqs. (4.1) and (4.4), we infer that the transition from a time-like to a space-like renormalized correlation function in position space amounts to the formal substitution ln it 2 → ln z 2 , (4.5) up to a possible change of the global sign.

Results
In this section, we collect our results for the case of standard QCD with the SU(3) gauge group and n f active quarks triplets. Notice that the results for the self-energy and the propagator of the "heavy" field, the ADs γ 2 and γ h as well as the corresponding RCs are gauge dependent. The expressions below are given in Landau gauge, as it is most relevant for lattice applications. Full results for a generic gauge group and including the gauge as well as the momentum/position dependence are appended in the arxiv submission of this paper as auxiliary files in computer readable format. Many results given in the text and in the auxiliary files have originally been obtained by other authors and have been included here for completeness. In particular the AD of the heavy-light current γ Q was computed at one, two and three loops in Refs. [54,55], Refs. [12,56] and Ref. [13], respectively. The AD γ h of the "heavy" field h v was computed at two and three loops in Ref. [12] and Refs. [13,57], respectively, and, recently, at four loops in Refs. [58,59]. The RCs Z 3 , Z 2 and Z a have been known through three loops for a long time [43][44][45]. The VEV Π in Eq. (1.1) was computed at two loops in Ref. [60] and at three loops in Ref. [61]. Notice that all these results for the ADs γ h , γ Q and for the VEV Π were obtained for the time-like choice of the vector v µ , with v 2 = 1.
Our contribution is to clarify the changes for the space-time choice v 2 = −1. We have also computed the VEV of the gluon off-light-cone operator in Eq. (1.2) in the threeloop approximation as well the corresponding anomalous dimensions at two loops. Our results are in agreement with the two-loop VEV found in Ref. [62] and the one-loop ADs first computed in Refs. [37,38].

Anomalous dimensions and renormalization constants
As follows from Eq. (3.4), ADs and RCs are not sensitive to the sign choice of v 2 . The analytic results for the ADs β, γ Q , γ h , γ 2 , γ ⊥⊥ and γ ⊥ are listed in Appendix A, and those for the RCs Z a , Z Q , Z h , Z 2 , Z ⊥⊥ and Z ⊥ in Appendix B.

Correlation functions: momentum space
Our result for the self-energy of the "heavy" field, Σ h , defined in Eq. (3.9) for the space-like case v 2 = −1 can be written as ¶ (δΣ s h ) Im n a n , where the first term corresponds to the time-like self-energy, 2) and the addenda δΣ s h arise when going over to the space-like case using the substitution rule in Eq. (3.43), Using Landau gauge, we find the following expressions: Notice that the heavy-field self-energy at v 2 = −1 acquires a large imaginary part, whereas the difference in the real part is minor. Next, we consider the momentum-space correlation function in Eq. (3.18). Our result forD reads (D t ) n a n + 2 n=0 (δD s ) Re n a n + i 2 n=0 (δD s ) Im n a n , (5.6) where, as above, the first term corresponds to the time-like correlation function, (D t ) n a n .

(5.7)
Notice that there is an overall minus sign between the time-like and space-like correlation functions, unlike for the self-energy Σ h .

Correlation functions: position space
As follows from Eq. (4.5), the time-like and space-like renormalized correlation functions in position space are given by identical expressions with the substitution it → z. Our result for the "heavy" propagator in position space is in agreement with Eq. (12) of Ref. [13] derived for v 2 = 1. For the correlation function in Eq. (1.1), we obtain where L z = ln (µe γ E z/2). For n f = 3, we get numerically: (F ) n f =3 = 1.0 + 7.05316 a s + 2.0 L z a s + 9.66546 a 2 s + 51.0988 L z a 2 s + 6.5 L 2 z a 2 s . (5.15) Notice that the higher-order coefficients in position space are considerably smaller than in momentum space, cf. Eq. (5.5).
Our results for the position-space functions Π ⊥⊥ and Π ⊥ , Π ⊥⊥ (z > 0, v 2 = −1) = 128 a z 4 1 + 2 n=1 (F ⊥⊥ ) n a n , Π ⊥ (z > 0, v 2 = −1) = − 128 a z 4 1 + 2 n=1 F ⊥ n a n , Nucleon matrix elements of these operators are usually called qPDFs and they are amenable to nonperturbative calculations on the lattice for space-like separations of the quark fields. At the same time, they are counterparts of similar time-like matrix elements that have been discussed in the past in the context of heavy-quark expansion in B-meson weak decays, and it is important to understand the relation between timelike and space-like renormalization and matrix elements.
We have shown, to all orders in perturbation theory, that the results for a generic Green function involving a qPDF operator at space-like and time-like separations are related by a specific substitution rule reflecting analytic continuation in the square v 2 of the four-vector v µ pointing along the Wilson line; see Eq. (3.43). The RCs and ADs are the same for space-like and time-like separations. This result is also relevant in the context of TMD factorization, where Wilson lines are shifted off the light cone to regularize rapidity divergences in TMD operators [16]. Our statement is that the ADs and RCs do not depend on the direction of the shift, space-like or time-like.
We have calculated the self-energy of the "heavy" field h v in the effective field theory of Eq. (2.1), the quark-antiquark qPDF AD and VEV, and all the RCs and ADs that are involved in the respective renormalization through three loops in the MS scheme. Our results agree with the literature as far as it goes. In addition, we have clarified the general RG pattern for the gluon qPDF operator in Eq. (1.2) and calculated its VEV through three loops, from which the two-loop ADs can be extracted avoiding pollution by gauge-noninvariant operators.
Our results can be used in lattice calculations aiming at the determination of quark and gluon PDFs, e.g., in the nucleon, if the linear UV divergences of lattice observables are removed by considering suitable ratios of matrix elements involving the same operator [34,35].