Quasi Transverse Momentum Dependent Distributions at Next-to-Next-to-Leading order

The partons’ transverse momentum can be explored with QCD lattice simulations by studying the quasi-transverse-momentum-dependent parton distribution functions (qTMDPDFs), which are factorized in terms of physical TMDPDFs and soft factors in the limit of the large hadron’s momentum. We present the next-to-next-to-leading order (NNLO) calculation of the coefficient function for this factorization. Together with already known expressions for anomalous dimensions, this result allows analysis of lattice data at NNLO perturbative accuracy

Introduction.Extraction of parton distribution and related quantities from the lattice simulations is a rapidly growing direction of QCD.The central tool for such extractions is the factorization theorems derived in the large hadron's momentum regime.Combining various equal-time operators and hadron states, one is accessing a variety of parton distributions [1,2].In this work, we study the so-called quasi-transverse momentum dependent (qTMD) parton distribution functions (PDFs) [2,3].The corresponding matrix element is defined as (we use the notation of Ref. [4]) where |P ⟩ is the (possibly polarized) hadron state with momentum P , q is the quark field, v and y are spacelike vectors (with y 0 = v 0 = 0), and b µ = y µ − v µ (vy)/v 2 .
The expression [a, b] represents a straight gauge link from point a to point b.The operator represents a quarkantiquark pair separated by y and connected by a stapleshaped Wilson line along direction v µ and of size L.
In the regime of the large momentum hadron P and the large length of gauge-link staple L, the qTMD matrix element (1) can be expressed via the physical TMD distribution through the factorization theorem, see Eq. (4), derived using various approaches in Refs.[4][5][6][7][8][9][10].A feature of the qTMD factorization theorem is that in addition to the physical transverse-momentum-dependent parton distribution functions (TMDPDF), it contains an extra function Ψ that accumulates the nonperturbative interaction between the parts of the staple gauge link (also called intrinsic soft factor [6]).The treatment of this function is slightly different in different approaches (compare, e.g., Refs.[7][8][9]), but conceptually the factorization theorem is the same in all cases.The main perturbative ingredient is the coefficient function, which is currently known at next-to-leading order (NLO) [8,11].
The operator (1) is localized in the equal-time plane, and thus the qTMD matrix element can be simulated within lattice QCD.The results of the simulations can be combined such that the functions Ψ cancel.In this way, one determines the nonperturbative Collins-Soper kernel, see Refs.[11][12][13][14][15][16].Alternatively, the Ψ function can be determined from the auxiliary procedure [6], and then one accesses the TMDPDF distribution [17].The precision of extraction crucially depends on the accuracy of the perturbative input, which is currently limited by the knowledge of the hard coefficient function.Importantly, this coefficient function is universal and is the same for all polarized quasi-TMDPDF of the leading power [8,9], and it is independent of the particularities of the nonperturbative definition of the internal soft factor.
Nowadays, the extractions of TMDPDFs from the data are routinely performed at next-to-next-to-leading order (NNLO) or N 3 LO order (see, for example, [18][19][20]), and recently were pushed to N 4 LO order [21].It has been demonstrated [18,20,22] that (at least) NNLO is required since the NLO is not sufficiently precise to describe the data from the modern experiments.Modern lattice simulations of quasi-transverse-momentumdependent parton distribution functions (qTMDPDFs) have yet to reach that order of precision.Still, nonetheless, the NNLO contribution is sizable since the typical scale of lattice simulations is about 1-3 GeV.In this paper, we present the expression for the hard coefficient function for the factorization of qTMDPDF at NNLO.Other perturbative ingredients of the factorization theorem (anomalous dimensions) are known at NNLO and higher.Therefore, using the result of this work, one could analyze the lattice data at complete NNLO.
Factorization theorem.The factorization theorem connecting qTMDPDFs and physical TMDPDFs is discussed in many articles [4][5][6][7][8][9][10], which we refer to for detailed discussion.Different spinor components of qTMD-PDF correlator (1) obey different kinds of factorization theorems [23].The projection to desired components is done by a contraction with appropriate Dirac matrix The most interesting cases are the components projected by Γ ∈ Γ + = {γ + , γ + γ 5 , iσ α+ γ 5 }, where n µ is a lightlike vector n 2 = 0 defined by the decomposition of hadron's momentum P µ = nµ P + + n µ M 2 /2P + (P 2 = M 2 ).For definiteness, we fix with v 2 = −1.The components projected by Γ ∈ Γ + obey the leading-power factorization theorem, where with and 11 is the physical TMDPDF of twist-two.The direction of the Wilson lines s is defined by the direction of the staple contour s = sign(L).In this way, different orientations of the staple contour give access to Drell-Yan or semiinclusive deep-inelastic scattering (SIDIS) definitions of TMDPDFs, which can be used to test their universality [24,25].We stress that the factorization limit is rather complicated (see the last term of (4)).Explicitly, it requires (vP ), L → ∞ at fixed-finite x and b.We also note that at this power accuracy, there is no difference between v − P + and (vP ) = P z , which is used as the hard scale.
The expression (4) is written in terms of renormalized functions.They are related to the bare functions as Here, the factor R renormalizes rapidity divergences (in most parts of schemes it is equal to the S −1/2 where S is the TMD soft factor).The factor Z J is the renormalization of the quark field in the axial gauge.The factors Z U 1 and Z Ψ1 are ultraviolet (UV) renormalizations of the (leading-twist) semicompact operators constituting the TMD distributions [23].Finally, the factor Z W depends on b and L and accumulates all divergent factors associated with the staple contour, such as power divergences of spacelike links [26], cusps divergences at point b + Lv and Lv, and other contributions [27].Importantly, the same divergences happen in the function Ψ, and thus we do deal with them in our computation of the coefficient function.
The structure of divergences cancellation in the qTMD factorization theorem is similar to those in factorization of Drell-Yan or SIDIS.So, the rapidity divergences cancel in-between soft factor contributions, Φ 11 and Ψ.The rapidity divergences leave no trace on the coefficient function, but introduce the rapidity scales ζ and ζ (see details Refs.[28][29][30]).The UV renormalization of TMD distributions cancels the infrared (IR) poles of the coefficient function.The cancellation happens only if The UV poles are renormalized by Z J 's and result in the overall scaling of the qTMDPDF operator.The scaling of functions ( 9) follows from their renormalization properties.For the Ω we have where γ J is the anomalous dimension of the heavy-light current, and γ W is the anomalous dimension associated with the renormalization of the staple link.The evolution of TMD distribution is [22,31] where Γ cusp is the anomalous dimension of the cusp of lightlike Wilson lines, and D is the Collins-Soper kernel.The Collins-Soper kernel is a nonperturbative function, which represents the interaction of light-quarks in the QCD vacuum environment [30].Finally, the Ψ function also evolves with the pair of equations where γ Ψ is the anomalous dimension associated with the finite part of the cusp anomalous dimension at the finite angle.The anomalous dimension 1 γ Ψ was computed at LO in Ref. 2 [8], and the NLO term was computed in this work.We found that γ Ψ coincides with the anomalous dimension associated with the heavy-quark [32] (v 2 > 0).This is not accidental, because the UV renormalization is insensitive to the sign of v 2 , as it is proven in Ref. [33].
The expressions for anomalous dimensions Γ cusp , γ J , γ V and γ Ψ are well known.At N 2 LO they can be found e.g. in Refs.[33][34][35].For the reader's convenience, we have collected all explicit expressions in Appendix A ( 29) - (32).The remaining anomalous dimension γ W and the Collins-Soper kernel are not important in the present work.
Quasi-TMD distribution.The qTMD distribution is an artificial construction that reduces to the physical TMD distribution in the asymptotic limit (vP ) → ∞ and L → ∞.Currently, there is no standard construction for this function, see discussion in Ref. [10].Probably, the most popular way [4,5,7,9,11] is to consider the function Using the factorization theorem (4), evolution equations and condition (10), one finds where dots indicate the power corrections explicitly given in the last term of Eq.( 4).Note, that in this formulation 1 By definition the anomalous dimension of the Ψ function is Evaluating it one should take into account that ζ ∼ µ 2 , because µ is the only dimensional scale of the Ψ function.Therefore, d ln(µ 2 /ζ)/d ln µ 2 = 0.For the detailed discussion see Ref. [4]. 2 Reference [8] provides an incorrect expression for LO γ Ψ .This mistake appeared due to the mismatch in definitions for the renormalization constant with earlier paper Z J ↔ Z −1 J .Once corrected, the expression for γ Ψ coincides with the one computed here or in Ref. [4].
the scaling equation for qTMD function is The nonperturbative part of the subtraction factor could be different in other constructions, and it does not affect C 11 .
Generally speaking, the scales µ and ζ are independent; therefore, one can define a more general function which reduces to (15) at ζ = µ 2 .This function satisfies the pair of equations Note that the evolution with respect to ζ has an opposite sign in comparison to ordinary TMD evolution (13).
Coefficient function.The qTMD operator can be written as a product J † v (y)ΓJ v (0), where the currents are Structurally, the current J v is similar to the renowned heavy-to-light current (see e.g.[36]), with v 2 = −1, and an open spinor index.The separation between currents y 2 ∼ b 2 is large in comparison to the hard scale (vP ) −1 , and thus any exchange of perturbative gluons between currents is power suppressed [4,[7][8][9].This essentially simplifies the problem of computation of C 11 and allows us to present it as where C 1 is the coefficient function for the factorization of a current (21) into the leading-twist semicompact operators.Because of this structure, the coefficient function is independent of the Γ once Γ ∈ Γ + .Therefore, it is universal for all eight leading-power components of the qTMDPDF matrix element.
Comparing the definitions (4), (9), and ( 22) we find that where we omit scaling arguments on the right-hand side for brevity.Note that the renormalization constant Z U 1 and C 1,bare are complex valued in such an approach [23].
The renormalization constants Z J , Z Ψ1 and Z U 1 are known at N 3 LO [33,35,37].For our NNLO computation, we took the expressions from the appendixes and auxiliary files of Ref. [33] (for Z J ) and Refs.[34,38] (for Z U 1 ), and reconstructed from known anomalous dimension [35] (for Z Ψ1 ).
The examples of diagrams contributing to C 1,bare are shown in Fig. 1.Note, that the same diagrams contribute to the computation of the matching coefficient of heavylight quark current in the heavy-quark effective theory (HQEFT) [39].The only difference between these computations is the sign of v 2 and that the momentum p is passing through the "light" quark line, while in the HQEFT matching coefficient computation, the momentum passes through the Wilson line.The computation is done in the dimensional regularization d = 4 − 2ϵ.The reduction to the base integrals is performed by the FIRE6 library [40].The result reads 1 + O(a 3 s ), (24) where and the expression for C (2) 1 is presented in Appendix A in Eq. (28).Combining C 1bare with the renormalization factors (and renormalizing a s ), we observe the exact cancellation of 1/ϵ poles.This provides a general check of the computation.
Substituting the renormalized expression for C 1 into Eq.(22), we obtain the NNLO coefficient function for the qTMD factorization theorem.It reads where L p is defined in Eq.( 8), C F = (N 2 c − 1)/2N c , C A = N c are eigenvalues of the Casimir operator of SU (N c ) algebra, N f is the number of active quarks, and ζ n is the Riemann ζ function.This expression is the main result of this paper.The NLO parts of the coefficient function and anomalous dimension γ Ψ coincide with the known results [8,11].
The logarithm part of the coefficient function can be derived from the evolution equations defined above.It satisfies Using explicit expressions for anomalous dimensions ( 29) -( 32) we confirm this.Note that using this equation and ( 26) one is able to compute the logarithm part of the N 3 LO coefficient function.We present it in Eq. (33).

Conclusion.
In this work, we have computed the coefficient functions for the factorization of the qTMD ma-trix at NNLO.We have checked the cancellation of poles between renormalization factors and coefficient function, which provides a check of the factorization theorem for the qTMD operator up to NNLO.These results also allow us to obtain the logarithm part of the N 3 LO coefficient function.The intermediate and final expressions are also attached to the publication in the Mathematica format.
In Fig. 2 we present the comparison of NLO, NNLO, and N 3 * LO coefficient functions at (vP ) = µ = 2GeV (i.e.L p = −2 ln x), where 3 * indicates that this coefficient function does not have the nonlogarithm term.At these energies, the coefficient function demonstrates a reasonable convergence for x > 0.2 (NNLO term provides ∼ 5% correction at most).Below x < 0.2 the convergence drops rapidly; e.g. at x = 0.1 the NNLO correction is ∼ 20%, and N 3 * LO is ∼ 40% (both in comparison to NLO).This shows the natural boundary x ≳ 0.2 for this approach.To access lower values of x one should find a way to improve the structure of the factorization theorem, either by resumming problematic logarithms (see, discussion on a similar problem for the pseudo-PDF case [41]) or by matching with a different type of factorization theorem).The coefficient function is just multiplicative, and thus it is straightforward to update the existing procedures including the NNLO correction.It is also important that the coefficient function is independent of the polarization quantum numbers.Therefore, all eight leading-power TMD distributions can be considered at the same NNLO precision.Let us mention that soon after the publication of the initial version of this manuscript the work [42] was published.In Ref. [42] authors derived the same coefficient function at NNLO by studying the threshold logarithms and related functions.Their results agree with ours.
Here, the expressions for Γ cusp and γ V are taken from Ref. [43], the expression for γ J is taken from Ref. [44], and the expression for γ Ψ is taken from Ref. [35].
Using these expressions and Eq. ( 26), together with Eq.( 27), we are able to determine the logarithmic part of the where c 3 is the unknown finite part.

Appendix B: Evaluation of diagrams
In this appendix, we provide extra notes about the computation of diagrams for the coefficient function at NNLO.The examples of diagrams are shown in Fig. 1.In total, there are 10 diagrams (including self-energy graphs).
The momentum enters the diagram via the quark line and goes out in the vertex.There is no momentum incoming into the Wilson line.For example, the second diagram shown in Fig. 1 reads , where d = 4 − 2ϵ is the parameter of dimensional regularization.It is important to keep track of +i0 prescriptions because incorrect prescriptions could lead to an improper sign of the resulting integral.For the leading-power computation, it is sufficient to consider p 2 = 0. Then the only dimensional parameter is ω = 2(P v).Also, only the good component (with respect to P ) of the quark field contributes to the leading-power term.To project the corresponding component we do where I is the diagram without a spinor multiplier.After projecting, the diagram decomposes into a sum of simpler scalar integrals, which have the general form F (a, b, c, d, e, f, g, h) = where all propagators have the +i0 pole prescription.Next, the integrals are reduced to the set of base integrals by the integration-by-parts relations (see Refs. [39,44] for a similar example).Specifically, we have used the FIRE6 library [40].Most parts of the base integrals are evaluated using successively one-loop integrals and are expressed with products of gamma functions.We found only two integrals with nontrivial topology.These integrals can be computed in the terms of higher-order hypergeometric functions or as an ϵ series [45] (for instance, we have used the Mellin-Barnes method).The results are F (0, 1, 0, 1, 1, 1, 0, 1) = (37) where γ E is the Euler-Mascheroni constant, and dots represent higher powers of ϵ series.Finally, by collecting the expressions together and expanding gamma functions in ϵ we obtain the bare expressions for each diagram.For example, the diagram (34) produces the following expansion:

FIG. 1 .
FIG. 1. Examples of diagrams contributing to C1 at NNLO.The momentum p enters through the quark line and exits from the quark-Wilson line vertex, as indicated on the left diagram.

FIG. 2 .
FIG.2.Comparison of NLO, NNLO, and N 3 * LO coefficient functions C11 as the function of x.The solid, dashed, and dotted lines represent the coefficient functions at NNLO, NLO, and N 3 * LO, correspondingly.The comparison is done for the value of (vP ) = µ = 2GeV, which is the typical setup for lattice computations.