Transverse Parton Distribution and Fragmentation Functions at NNLO: the Gluon Case

We calculate in this paper the perturbative gluon transverse momentum dependent parton distribution functions (TMDPDFs) and fragmentation functions (TMDFFs) using the exponential regulator for rapidity divergences. We obtain results for both unpolarized and linearly polarized distributions through next-to-next-to leading order in strong coupling constant, and through ${\cal O}(\epsilon^2)$ in dimensional regulator (finding discrepancy for the linearly polarized gluon TMDPDFs with a previous result in the literature). We find a nontrivial momentum conservation sum rule for the linearly polarized component for both TMDPDFs and TMDFFs in the ${\cal N}=1$ super-Yang-Mills theory. The TMDFFs are used to calculate the two-loop gluon jet function for the energy-energy correlator in Higgs gluonic decay in the back-to-back limit.

In this paper, we present the results for the perturbative matching coefficients at Next-to-Next-to Leading Orders (NNLO) for gluon TMDPDFs and TMDFFs. They are relevant to the production and decay of the Higgs boson, top quark pair, and J/ψ at low transverse momentum at NNLO or Nex-to-Next-to-Next-to Leading Logarithms (N 3 LL) order. We also provide the bare results at NNLO through O( 2 ), which do not contribute to the renormalized TMD coefficients at this order, but are relevant for future N 3 LO calculation. It is well-known that, direct calculation of the TMD matching coefficients requires some form of regularization in addition to the usual dimensional regularization [1,5,24,[37][38][39][40][41][42][43][44]. In this paper, we adopt the exponential regularization scheme for the rapidity divergences [42].
For the gluon TMD coefficients, the non-trivial spin structure leads to two independent tensor structures in transverse impact parameter space. They are known as the unpolarized coefficients and linearly polarized coefficients [45]. The linearly polarized coefficients arise from helicity-flip contribution in gluon splitting, and are suppressed by one power of α s compared with the unpolarized ones. Their contributions to physical cross section have been discussed for diphoton production [46], Higgs production [15,[47][48][49], quarkonium production [50][51][52], γ * plus jet production [53,54], heavy quark pair [55] and dijet production [56]. In this paper we present results through O(α 2 s ) (two loops) for both polarizations. For the linear polarization, this is formally at the NLO accuracy, since its LO contribution already starts at O(α s ). We note that with the two-loop linearly polarized contributions presented in this paper, the linearly polarized ingredients needed for N 3 LO calculation [19,57] for Higgs production using the q T subtraction formalism [58] are completed. The reason is that the pure three-loop linearly polarized component will not contribute to the cross section at O(α 3 s ), since the tree-level linear polarized contribution vanishes, and its interference with the unpolarized contribution vanishes.
As a by-product of this calculation, we find an interesting momentum conservation sum rule for the linearly polarized TMD coefficients in the N = 1 supersymmetric limit. The sum rule imposes non-trivial constraint to the gluon-to-gluon and quark-to-gluon TMD coefficients, and are found to be satisfied by our two-loop results.
The unpolarized gluon TMDPDFs have been computed before using different regulators for the rapidity divergences. Our results can be compared with them when combining with the rapidity-regulator dependent TMD soft function properly [59]. We have done this exercise and found complete agreement with Refs. [60,61]. For unpolarized gluon TMDFFs, results have also been given through NNLO [62]. Our results agree with them for most of the terms, except for a term of the form C 2 A π 4 δ(1 − z) in the gluon-to-gluon TMD coefficient. We have also found similar discrepancy in the quark-to-quark coefficient as reported in Ref. [63]. Very recently, the two-loop results for the linearly polarized gluon TMDPDFs have been given in Ref. [64]. We have compared our results with Ref. [64] by constructing rapidity-regulator independent TMDPDFs, and found discrepancy with Ref. [64]. We have performed several checks to our calculations, which will be explained below.
This paper is organized as follows. In Sec. 2, we give the bare and renormalized results for the gluon TMDPDFs through NNLO. In Sec. 3, we give the bare and renormalized results for the gluon TMDFFs through NNLO. In Sec. 4 we use the gluon TMDFFs to calculate the gluon jet function, which is relevant for the Energy-Energy Correlations for Higgs gluonic decay in the back-to-back limit. We conclude in Sec. 5. We collect the relevant perturbative ingredients in the appedix A.

Gluon TMDPDFs
The bare gluon TMDPDF can be defined in terms of SCET [65][66][67][68][69] collinear gauge fields where A a,µ n⊥ is the gauge invariant collinear gluon field with color index a and Lorentz index µ. For sufficiently small b ⊥ , the gluon TMDPDFs admit operator production expansion onto the usual collinear PDFs, where the summation is over all parton flavors i. The coefficient functions can be decomposed into two independent Lorentz structures, where we have defined two scalar form factor I bare gi and I bare gi , which can be projected out using 2) do not depend on the actual hadron N . In actual calculation, one can replace the hadron N with a partonic state j. Furthermore, the usual bare partonic collinear PDFs are just φ bare i/j (x) = δ ij δ(1 − x), so one has (2.5) The TMDPDFs, as well as their matching coefficients, contain both UV and rapidity divergences. We adopt dimensional regularization for the UV, and exponentional regularization [42] for the rapidity divergences. In the following subsection, we present the bare results for the coefficient functions through two loops in QCD.

The bare results
The relevant diagrams for gluon TMDPDF through two loops are generated with the code Qgraf [70]. We use an in-house Mathematica code to substitute in the SCET Feynman rules, and use Form [71] to carry out necessary color and algebra manipulation. We employ reverse unitarity [72] to convert phase space integral into loop integral for the purpose of integral reduction. We use Fire5 [73] and LiteRed [74] to reduce the integrand into the so-called master integrals by Integration-By-Parts identities [75]. The resulting master integrals have been solved in Ref. [63] using differential equation method [76,77]. For the details we refer to Ref. [63]. The bare coefficient functions can be expanded in terms of QCD bare coupling as and similarly for the form factor I bare gi and I bare gi , where S = 4πe −γ E µ 2 0 /µ 2 . The oneloop bare coefficient functions up to O( 2 ) are where we have defined and with b 0 = 2 e −γ E . Our notations are the same as in Ref. [63]. The two-loop bare results up to O( 0 ) are where all end point divergences ln k (1 − x)/(1 − x) should be taken as plus-distribution, which is defined as for a smooth test function g(x). We write our results in terms of Harmonic PolyLogarithms [78]. We use the Mathematica package Hpl [79] to manipulate them. We use the standard shorthand notation H a 1 ,...,an ≡ HPL(a 1 , . . . , a n ; x) . (2.12) We have also obtained the two-loop results through O( 2 ), which are required for future N 3 LO calculation. We do not show them here to save space, but instead include them as ancillary file in the arXiv submission of this paper.

Renormalization counter terms and zero-bin subtraction
The bare results contain both UV and rapidity divergences to be renormalized, while the rapidity divergence is already renormalized within exponentional regularization upon the substitution τ → 1/ν. To proceed, we first perform the usual coupling constant renormalization in MS scheme, where γ E = 0.5772 . . . is the Euler constant, µ 0 is the mass parameter in dimensional regularization, α s = α s (µ) is the renormalized coupling evaluated at the renormalization scale µ. The one-loop QCD beta function for N f light flavor is given by Next, we perform a zero-bin subtraction [80] to remove the overlapping contributions between collinear and soft modes. For the purpose of computing perturbative matching coefficients, we can used the dimensional regularized collinear PDFs in a partonic state, where P ij (x) are space-like splitting kernel, whose explicit expressions are collected in the appendix. After this, the remaining UV divergences in the TMDPDFs can be removed by a multiplicative renormalization counter term. These steps can be summarized as where B bare,µν g/j and S 0b (α s ) are the bare TMDPDFs and bare zero-bin soft function in terms of renormalized coupling α s , Z B g is the multiplicative operator renormalization constant, and in the last equality we use the renormalized version of operator product expansion in Eq. (2.2). The zero-bin soft function can be found in appendix. A.6.

Renormalized coefficient functions
In this subsection, we present the detailed results for the renormalized matching coefficient functions through O(α 2 s ). The renormalized coefficients obey a RG equation and a rapidity evolution equation [39] (2.18) The relevant anomalous dimensions are collected in Appendix. A.1. The tensor decomposition for the renormalized quantities is the same as given in Eq. (2.3). Throughout this paper, we define the pertubative expansion according to Then, the renormalized scalar form factors are given by gi (x) + I gi (x) , (2.20) where we have shown explicitly the scale-dependent part and scale-independent part. At each order, the scale-dependent part are determined by RG equations and universal anomalous dimensions. They serve as strong check of the results from Feynman diagram calculation. The genuine new results of direct calculations are the scale-independent terms. At one loop they are given by where r b (x) is defined in Eq. (2.8). At two loops, we find The two-loop scalar form factors corresponding to the first tensor structure, I gi , have been computed for a while [60,62]. Results for the second tensor structure I (2) gi , also known as the linearly polarized contribution, appeared very recently, using a different rapidity regulator [64]. Although the calculations are performed with different rapidity regulators, the results can be compared by constructing a rapidity-divergence free combination of TMDPDFs. In our case, this can be done by multiplying the renormalized coefficient functions with the square root of the TMD soft function, for i = q, g, and then expand in α s . The renormalized TMD soft function can be found in Ref. [59]. The explicit two-loop results for the renormalized soft function are collected in Eq. (A.14) of appendix. A.5. We find that for the first tensor structure, our two-loop results are in full agreement with those in the literature [62,81]. However, we find substantial difference for the two-loop linearly polarized results with those presented in Ref. [64]. We note that the two-loop results in Ref. [64] contain transcendental functions and zeta value up to weight 3, namely Li 3 and ζ 3 , while in our results, I gg and I (2) gq in Eq. (2.22), only weight 2 functions and zeta values are presented. Since the linearly polarized contribution first appears at one loop with rational functions only, our results are in agreement with the expectation that transcendental weight only increases by 2 at each loop order. In the next subsection, we present a N = 1 supersymmetry sum-rule for the linearly polarized gluon contribution, which provides further check to our results.

N = 1 supersymmetry sum rule for the linearly polarized gluon contribution
Our explicit two-loop calculation reveals an interesting momentum conservation sum rule for the linearly polarized gluon contribution in the N = 1 supersymmetric limit. This limit is obtained from our results by setting C A = C F = N f and using T F = 1/2. The sum rule is then written as At one-loop, the sum rule is satisfied trivially, since I (1) , see Eq. (2.21). At two loops, we find substituting this into Eq. (2.24) we indeed get zero. While it is well-known that the splitting functions obey momentum conservation sum rule, in general the sum rule breaks down for the matching coefficient functions, which are cross-section level quantities. Therefore, the sum rule in Eq. (2.24) is somewhat surprising. It would be interesting to see if it continues to hold at three loops. It would also be interesting to have a structural understanding of it. Since the sum rule is nontrivial, it also provides strong check to the two-loop results for the linearly polarized gluon contribution presented in this paper.

Gluon TMDFFs
To specify the definition for the gluon TMDFFs, it is necessary to specify a reference frame first. In the hadron frame, where the detected hadron has zero transverse momentum, the gluon TMDFFs can be defined as where N (P ) has zero transverse momentum. In actual calculation, it's also convenient to define the fragmentation functions in the parton frame, where the parton which initiates the fragmentation has zero transverse momentum. The parton frame TMDFFs are related to the hadron frame ones by where we denote the bare TMDFFs in the parton frame by F bare,µν The parton frame has the advantage that the renormalization counter terms are slightly simpler compared with the hadron frame, i.e.

(3.3)
We refer to Refs. [1,63] for more detailed discussion on the difference of the reference frames.
To compute the gluon TMDFFs, we can use the methods as in the calculation of TMDPDFs. Alternatively, we can exploit crossing symmetry to obtain the results from TMDPDFs. We have performed calculation in both ways, and find the same results. In the next subsection, we give some details on the calculation based on crossing relation.

Bare gluon TMDFFs from crossing
To explore the crossing symmetry between the TMDPDFs and TMDFFs, we write down their definitions in momentum space. For simplicity we first consider the g → g case, for TMDPDFs. In Eq. (3.4), |M B µν gg (p 1 , l j , k i )| 2 is the squared amplitudes for space-like g → g splitting, with multiple real radiations (k i ) or virtual momentum exchange (l j ), and p 1 is the momenta of the initial-state gluon entering hard scattering. For TMDFFs in hadron frame, we have where |M F µν gg (p 2 , l j , k i )| 2 is the squared amplitudes for time-like g → g splitting, and p 2 is the momenta of the final-state detected gluon, which has zero transverse momentum in hadron frame. The squared amplitudes in the integrand for TMDPDFs and TMDFFs are related through the following crossing relation, It is not difficult to see that that is, where we have defined The rule Q 1 → −Q 2 is important for resolving the ambiguity from analytic continuation of ln(1 − x) (see also discussions in [82][83][84][85][86] for analytic continuation of splitting funtions), which can have two possibilities, We argue that one should set κ = 0 for ln(1 − x) terms originate from rapidity divergences, and κ = 1 for ln(1−x) terms originate from virtual corrections. For virtual corrections, the analytical continuation is unambiguous, since it is determined by Feynman's iε prescription.
To understand the prescription for the ln(1−x) from rapidity divergences, we first introduce the following dimensionless variables, Using these variables, the exponential regulator in TMDPDFs and TMDFFs can be separately written as [42,63] where we have used that 2k 0 =n · k + n · k. These two expressions can be exactly related to each other through the following rules, The exponential regulator regularizes the rapidity divergences in the phase space integral only. Therefore, the ln(1 − x) terms that originate from rapidity divergence should not develop imaginary part under the replacements x → 1/z, Q 1 → −Q 2 . This is indeed the case since the logarithm always appear in a specific combination (3.14) In practical calculation, one prescription to resolve the ambiguity in Eq. (3.10) is to use the the following crossing rules (3.15) Moreover, taking into account the relation between parton frame and hadron frame, which is given in Eq. (3.2), and the color and spin factors, the analytic continuations from TMD-PDFs to TMDFFs reads In the final step, we take the real part in Eq. (3.16). Note that the analytic continuations in Eq. (3.16) are valid in the region The contributions from the end point, δ(1 − x) and δ(1 − z), are invariant under crossing. The analytic continuations in Eq. (3.16) can be also used to extract time-like splitting functions from space-like ones. It is similar to the analytic continuations of splitting function performed in Refs. [84,85]. There the calculations start from the crossing of unrenormalized structure functions and work well for two-loop splitting functions. At three loops, the direct analytic continuations cause some issues (see [84,85] for details), so it is expected our procedure presented here may also cause similar issues at three loops. However, this is beyond the scope of this paper.

Renormalization counter terms and zero-bin subtraction
The renormalization of the TMDFFs are similar to the renormalization of the TMDPDFs. In the parton frame, we use the following dimensional regularized collinear FFs for the counter terms: where P T ij (z) are the time-like splitting kernel, whose explicit expression through two loops can be found in appendix. A.3. In hadron frame, the counter terms from collinear FFs involve additional factor of z −2 [62], which arises from the phase space factor. The zerobin subtraction also follows closely the TMDPDFs.

Renormalized coefficient functions
In this subsection, we present the renormalized coefficient functions through two loops for both polarizations. We separate the results into scale-independent part and scaledependent part. The scale-dependent part is determined by the the RG equation (3.20) and rapidity evolution equation (3.21) The two polarization form factors can be extracted through Compared with Eq. (2.9), for TMDFFs we use a slightly different definition for L Q , The explicit solutions to both equations through to two loops are given by ig (z) .

(3.24)
Our explicit diagrammatic calculation or analytical continuation reproduces all the scaledependent terms, which serve as a first check to the calculation. The genuine new results at each order are the scale independent part. At one loop they are given by where we have defined The two-loop results are qq . Like for quark TMDFFs, we have computed the gluon jet function for the EEC in the back-to-back limit and tested with known constraints to verify our results. We shall present the details of this check in the next subsection.
Our two-loop results for the linearly polarized coefficients are new. Similar to the TMDPDFs, the two-loop linearly polarized contributions involve only weight 2 functions and zeta values. Furthermore, they obey a N = 1 supersymmetric momentum-conservation sum rule, (3.28) We note that there is a change of sign in the sum rule for C qg compared with the corresponding sum rule for TMDPDFs in Eq. (2.24).

Gluon jet funcion for the EEC in the back-to-back limit
The gluon TMDFFs can be used to calculate the gluon jet function for the back-to-back resummation of the EEC in the Higgs gluonic decay [87]. Such jet functions also appear in the back-to-back resummation for TEEC [36]. The gluon jet function can be expanded as The expansion coefficients through two loops are given by the second moment of the gluon TMDFFs, where the two-loop scale-independent constant for gluon jet is The gluon jet constant in Eq. (4.3) enters the back-to-back contact terms in the Higgs EEC at two loops. These constants provide important ingredients for N 3 LL resummation of the EEC in Higgs gluonic decay in the back-to-back limit. They have also been used in Ref. [88] to extract the collinear contact terms for the Higgs EEC using energy-conservation sum rule.
We note that the two-loop jet function has also been computed explicitly [89] from taking the asymptotic expansion of the differential equation satisfied by the master integrals for the Higgs EEC at finite angle [87]. Our results in this paper for the jet function agree perfectly with the brute force calculation in Ref. [89]. Another independent check of our results comes from the momentum-conservation sum rule for the EEC, discovered in Refs. [90,91]. A similar discussion for the quark jet function can be found in Ref. [63].
We stress that the jet functions are the second moment of the sum of all the TMDFFs, including both the unpolarized and linearly polarized contributions. Therefore, the aforementioned checks also apply to all the TMDFFs presented in this paper.
It is also interesting to consider the contribution to gluon jet function from linearly polarized gluon alone. They are given by J g,l.p.
It's interesting to note that the linearly polarized gluon contribution in Eq. (4.4) vanishes in the N = 1 supersymmetry limit, where C F = C A = N f . This is simply the consequence of the supersymmetric momentum-conservation sum rule in Eq. (3.28).
Recently, EEC in the back-to-back limit has also been studied at hadron collider, the so-called TEEC. Due to the special geometry at hadron colliders, the linearly polarized gluon does not contribute to TEEC [36]. In that case, the gluon jet function is simply given by Therefore, our results also provide the missing ingredients for the N 3 LL resummation of the TEEC in the back-to-back limit.

Conclusion
In this paper we have presented the perturbative gluon TMD coefficient functions through two loops, for both unpolarized and linearly polarized coefficients. Our calculation was performed using the exponential regulator for the regularization of rapidity divergences. We have obtained results through O( 2 ) at two loops, which are relevant for future three-loop calculation. Our results are rapidity regulator dependent. Rapidity regulator independent results can be obtained by multiplying the TMD soft function properly, see Eq. (2.23). We have compared our results with those in the literatures, and found perfect agreement in most cases. However, discrepancy with Ref. [64] was found for the linearly polarized gluon TMDPDFs at two loops, and with Ref. [62] for the δ(1 − z) term in the gluon-to-gluon TMDFF. As a by-product of this calculation, we found a momentum conservation sum rule in the N = 1 supersymmetric limit for the linear polarization of gluon TMDPDFs and TMDFFs. Our results provide important ingredients for the precision studies of gluon TMDs in collider experiments.

A Anomalous dimensions, splitting functions, renormalization factors and the TMD soft function
In this Appendix, we list some necessary ingredients which enter our calculation.

A.1 Anomalous dimensions
For all the anomalous dimensions entering the RGEs of various TMD functions, we define the perturbative expansion according to where the coefficients up to O(α 2 s ) are The cusp anomalous dimension Γ cusp can be found in [92,93] and rapidity anomalous dimension γ R can be found in Refs. [59,94]. The hard and soft anomalous dimensions γ H and γ S can be extracted from the two-loop gluon form factor [95], and can also be found in, e.g., Refs. [96,97]. Finally, the beam anomalous dimension γ B is related to γ S and γ H through γ B = γ S − γ H . The QCD beta function is defined by with [98][99][100][101][102] A formula particularly useful for us is

A.5 Renormalized TMD soft function
The exponentially regularized TMD soft function is given by [59] S gg (b ⊥ , µ, ν) = exp α s 4π c s where the scale-dependent terms are 15) where L R = L ⊥ + L ν with L ν = ln(ν 2 /µ 2 ). The scale-independent terms are The bare zero-bin soft function is equivalent to bare TMD soft function, where we need to expand c s 1 and c ⊥ 1 to order O( 2 ),