Linearly polarized gluons at next-to-next-to leading order and the Higgs transverse momentum distribution

We calculate the small-b (or large-qT) matching of transverse momentum de- pendent (TMD) distribution for linearly polarized gluons to the integrated gluon distributions at the next-to-next-to-leading order (NNLO). This is the last missing part for the complete NNLO prediction of the Higgs spectrum within TMD factorization. We discuss the numerical impact of the correction so derived to the qT -differential cross-section for Higgs boson production and to the positivity bound for linearly polarized gluon transverse momentum distribution.


Introduction
The gluon-gluon fusion is the leading channel for the Higgs boson production in hadronhadron collisions [1][2][3]. The transverse momentum dependent (TMD) factorization of Higgs production has been demonstrated to follow the same pattern as the Drell-Yan/vector boson case in different frameworks [4][5][6][7][8] and in this sense it has been reviewed in [9]. Within the TMD factorization theorem, which describes the Higgs production at small transverse momentum, there are two dominant terms in the factorized cross-section. Those terms correspond to the fusion of unpolarized and the linearly polarized gluons [10][11][12]. Schematically, it reads where σ gg→H is the factorized gluon-gluon-Higgs cross-section, x A,B are the collinear fractions of gluon momenta, f 1 is the unpolarized gluon traverse momentum dependent parton distribution function (TMDPDF) and h ⊥ 1,g is the linearly polarized gluon TMDPDF (lpT-MDPDF) that was proposed as an independent distribution a long ago by Mulders and Rodrigues [13].

JHEP11(2019)121
In TMD factorization each TMD distribution (f 1,g and h ⊥ 1,g in this case) is an independent fundamental non-perturbative function. In order to sensibly construct a TMD for any practical purpose it is fundamental to include the asymptotical small-b limit, where each TMD distribution match to collinear parton distributions and the matching coefficient is calculable in QCD perturbation theory [5,6,14]. The modern state-of-the-art of perturbative calculations these matchings is the next-to-next-to-leading order (NNLO) of perturbative series, see [15][16][17][18]. Such a high order is required because of the sizes of theoretical and experimental uncertainties, see e.g., [19,20]. Also, it is required for the use of NNLO TMD evolution, which is necessary to perform an accurate global analysis of highand low-energy data [21,22]. The small-b limit of the unpolarized gluon TMDPDF, f 1 , has been calculated at NNLO in [15,16]. However the small-b limit of the lpTMDPDF, h ⊥ 1,g is known only at one-loop [9,12,17] and as such it has been used in ref. [23]. 1 In this work, we fill this gap, providing the calculation of h ⊥ 1,g at two loops and estimating the impact of this correction on the Higgs transverse momentum spectrum. The calculation can be performed using the same techniques as in refs. [16,[24][25][26].
The necessity of the present calculation comes from the fact that the perturbative counting in TMD formalism is slightly different from the one used in resummation approach. In fact, using standard resummation, see e.g. [12,20,[27][28][29][30], the small-b expansion is incorporated into the factorization formula, ignoring the non-perturbative TMD effects and one worries only about the perturbative expansion of the cross section. The TMD factorization includes the resummation for large enough q T , however one has different requirements in the realization of the perturbative series. So, while in the usual resummation the whole bracketed factor in eq. (1.1) should be given at a certain perturbative order, in TMD factorization each distribution should be matched independently to its collinear counterpart at the same given order. Both approaches are consistent with computing the small-b expansion at the same order. The case of linearly polarized gluon contribution to eq. (1.1) is special because the tree-level matching accidentally vanishes. The counting of perturbative orders in TMD factorization is reported later in the text.
The result obtained in this work is relevant for many cases beyond the Higgs boson production. In particular, there are processes that are also sensitive to lpTMDPDF and that are addressed in the literature [31][32][33][34][35]. Among these it is worth a special mentioning the case of heavy-quark production [36][37][38][39][40][41][42], which is relevant at LHC, future Electron-Ion Collider (EIC) or the LHeC. Another important topic is the positivity bound for gluon TMDPDF derived in [13], This positivity bound is expected to saturate at small-x due to the McLerran-Venugopalan model [43]. Our calculation shows that this bound is easily violated by loop corrections but could be restored by non-perturbative corrections. In this way, the relation in eq. (1.2) could be considered as a strong restriction on transverse momentum dependence of partons.

JHEP11(2019)121
The two-loop calculation presented here is structured in a way similar to the case of unpolarized gluons, evaluated in [16]. We find it sufficient to recall the basic principles and notation in section 2, which can be skipped by the reader already acquainted with topical works. The computation has requested the calculations of several new master integrals which are reported in the appendix. The final result for the NNLO matching of h ⊥ 1,g onto collinear gluon PDF is presented in section 3. The NNLO matching calculated here has been incorporated into artemide [44, 45], which was used to perform a qualitative numerical estimation of lpTMDPDF to Higgs-production cross-section at NNLO-N 3 LL. The results of the phenomenological analysis are discussed in section 4.

Definition
The TMD distribution of gluons in a hadron is given by the following matrix element where n is a light-like vector, F µν is the gluon field strength tensor, andW denotes the half-infinite Wilson line in the direction ñ The Wilson linesW n are taken in the adjoint representation of the gauge group. We use the standard notation for the light-cone components of vector v µ = n µ v − +n µ v + + g µν T v ν (with n 2 =n 2 = 0, n ·n = 1, and g µν T = g µν − n µnν −n µ n ν ). The decomposition of the gluon TMD distribution over independent Lorenz structures contains 8 components [9,13]. Two of these structures survive in the case of unpolarized hadron For future necessity, the decomposition in eq. (2.3) is given in d = 4 − 2 -dimensions as it was defined in [15,17]. Both f 1 and h ⊥ 1 contribute to the gluoninduced TMD processes on equal foot. Although these functions share some common properties, they are completely independent non-perturbative functions that are to be extracted from the experiment.
The usage of a d-dimensional definition for the decomposition in eq. (2.3) is important for the following two-loop calculation because the -dependent parts influence the result. The definition in eq. (2.3) is the standard one [15,17] written such that the unpolarized part coincides with the standard definition of the unpolarized TMDPDF, see e.g. [5,15,16]

JHEP11(2019)121
(here dots denote the staple gauge link, as in (2.1)), (2.4) whereas the linearly-polarized tensor is orthogonal to it. In turn the lpTMDPDF is given by Sometimes, one would like to use TMD distributions defined in the momentum space. The relation between coordinate and momentum representation is the usual one [9,13] (here in d = 4 dimensions),

OPE at small-b
At small-b the TMD operator can be matched to the collinear operators by means of operator product expansion (OPE). This relation is important because it constrains the model for TMD distributions at small values of b. Moreover, at large values of Q, where the TMD evolution factor significantly suppress the large-b part of the Fourier integral, the small-b OPE provides the dominating input to the cross-section (see e.g. [9,12,23] for studies related to Higgs boson processes). The systematic description of the small-b OPE applied to TMD operators can be found in ref. [46]. In the present case, it results into the following expressions where the sum runs over the active parton flavors (quarks and gluon), and f 1 (x, µ) is unpolarized collinear distributions defined as usual

JHEP11(2019)121
Concerning the notation, here and in the following we distinguish the unpolarized TMD-PDF f 1 (x, b) and unpolarized collinear PDF f 1 (x) by the number of arguments. The scales µ and ζ in eq. (2.9), (2.10) are the scales of TMD evolution discussed in the next section. The scaleμ is the scale of OPE, that is not related to the TMD evolution scales and whose dependence cancels in the convolution of coefficient function and collinear distribution. The coefficient functions (also known as matching coefficients [5]), C and δ L C, are to be calculated in QCD perturbation theory. The three-order calculation yields where a s = g 2 /(4π) 2 is QCD coupling constant. Nowadays, the coefficients C f ←h (x, b) are known at a 2 s -order (NNLO) [15,16,24,47], whereas coefficients δ L C f ←h (x, b) are known at a s -order (NLO 2 ) [9,12,17]. In the following section we present NNLO expression for δ L C g←f , which allows to consider these distributions at the same level of accuracy.
The corrections to the OPE at higher powers of b are unknown but at large value of b 2 the OPE becomes divergent. Thus, in practice, for the description of the TMD distributions one typically uses a phenomenological ansatz that matches the OPE results at small-b to a non-perturbative input at large-b. It can be written in the form and a similar expression can be used for f 1 (x, b) with a different f 1NP (x, y, b 2 ) and the corresponding matching coefficient. In eq. (2.15) we omit scale variables, and the function h ⊥ 1NP is an arbitrary function with the only constraint which is necessary to be consistent with the small-b limit of the TMD. A similar ansatz has been used also for the quark TMD, and the respective non-perturbative correction has been called f N P in refs. [21,22]. Up to now the non-perturbative correction to the quark TMD is the only one which has been extracted from data. In order to have some phenomenological result here we also choose We comment about the consistency of this choice in section 4.

Renormalization of TMDPDF
The TMD operator contains ultraviolet (UV) and rapidity divergences. Both these divergences can be renormalized (the all-order proof of renormalization for rapidity divergences is given in ref. [8]) by the corresponding renormalization factors. Hence, the renormalized (or physical) TMD distribution depends on two scales µ (the UV renormalization scale)

JHEP11(2019)121
and ζ (the rapidity divergences renormalization scale). The renormalized expression for the TMD distribution Φ g←h reads where Φ µν;unsub. g←h denotes the bare or unsubtracted TMD distribution, either f 1 either h ⊥ 1 , since the TMD renormalization is independent of polarization properties. In eq. (2.17) we present explicitly the dependence on regularization parameters: is the parameter of dimensional regularization (d = 4−2 ) that regularizes UV divergences, that are renormalized by the factor Z g ; δ is the parameter of δ-regularization [16,25] which regularizes rapidity divergences that are renormalized by the factor R g . The renormalization factors Z g and R g are ordered such that the renormalization of rapidity divergences is made before to the renormalization of UV divergences as it was done in similar NNLO calculations [16,24,26]. The final result is independent of the subtraction order.
The rapidity renormalization factor can be related to the TMD soft factor [8], which is the vacuum expectation value of certain Wilson loop [5,6,8], whereW n andWn are Wilson lines along n andn (2.2). In the case of gluon operators the Wilson loop is in the adjoint representation. The rapidity divergences are regularized by the δ-regularization, which consists in suppression of the gluon field in a Wilson line by exponential factor, A + (nσ + x) → A + (nσ + x)e −δ|σ| . The rapidity divergences reveals as ln(δ). In this scheme the rapidity renormalization factor is [8,25,48] The variable p + is parton momentum [46], and is required to define the Lorentz invariant scale ζ. Note, that the definition (2.19) also contains finite at δ → 0 terms, which can be seen as a scheme-dependence. Commonly, the scheme dependence is fixed by condition that no remnants of the soft factor appear in the hard part of the factorization theorem [5,8]. Definition (2.19) satisfies this condition. The UV renormalization factor is taken in MS-scheme. The (µ, ζ)-dependence of gluon TMD distribution is provided by a pair of evolution equations known up to three-loop order inclusively [49][50][51][52]. Note, that the renormalization factor Z g also contains the gluon-field renormalization part, therefore, where the symbol AD extracts the coefficient of −1 with a pre-factor n! at the n th perturbative order. Anomalous dimensions γ g and D satisfy the integrability condition (also known as Collins-Soper equation [53]) where Γ cusp is anomalous dimension for cusp of two light-like Wilson lines (in the adjoint representation). Due to this equation the expression for γ g can be rewritten in the form where γ g V is anomalous dimension of the vector form factor for gluon. The rapidity anomalous dimension D g has not such a simple representation due to the presence of an extra dimensional parameter b 2 . It generally contains all powers of logarithms ln(µ 2 b 2 ), that at some large values of b 2 turns to some non-perturbative function [54].
Due to the integrability condition in eq. (2.23) the system of evolution equations in eq. (2.20), (2.21) has a unique solution: where the TMD renormalization factor reads Here, P is arbitrary path in (µ, ζ)-plane connecting (µ 1 , ζ 1 ) and (µ 2 , ζ 2 ). The eq. (2.26) is in principle independent of the path P , however the truncation of the perturbative series makes some choices more preferable, for the detailed discussion see ref. [19]. In particular, in section 4 we use the special practically-convenient path that corresponds to ζ-prescription introduced in [19,21]. We again stress that the TMD evolution equations and their solution of eq. (2.25) do not depend on the polarization, and thus it is exactly same for unpolarized TMDPDF f 1 and lpTMDPDF h ⊥ 1 .
3 Matching coefficient for lpTMDPDF at NNLO

Evaluation of the matching coefficient
The coefficient function for OPE at twist-2 level can be deduced from the calculation of matrix elements with free parton states with subsequent matching of the result on the desired OPE structures eq. (2.10). Therefore, the task is naturally split into two steps:

JHEP11(2019)121
the evaluation of parton-matrix element and the matching. This procedure is well-known, see e.g. [5,6,9,16,24,26], in this section we present only minimal details and specifics of calculation of lpTMDPDF. The evaluation of parton matrix elements of the TMD operators at two-loop level is the most complicated part of the present work. We have used the same technique that was used by our group for NNLO evaluations in refs. [9,16,26], where we refer for extra details. In the case of lpTMDPDF the main complication comes from the rich vector structure, which is reduced to scalar products by projection factor in eq. (2.5), and the use of unpolarized parton states with momentum p µ = p + n µ . In this aspect the current computation is similar to evaluation of the pretzelosity distribution [26] albeit with significantly larger number of loop-integrals. The reduction of integrals to master integrals and some details of their evaluation is presented in the appendix A.
The outcome of each diagram at NNLO has a generic form The functions g 2 and g 3 exactly cancel in the sum of all the diagrams (and this fact can be also traced in the sum of sub-classes of diagrams) because they represent IR divergences. The last two terms represent the rapidity diverging pieces, and thus the functions g 4 and g 5 are canceled by the rapidity renormalization factor. However, due to the absence of three-order term, the functions g 5 cancel in the diagrams. The cancellation of all these pieces provides a check of the calculation. Summing together the diagrams we obtain the un-subtracted expression for TMDPDF on free-gluon states. Let us introduce the notation for perturbative series where a s = g 2 /(4π) 2 . The tree-order term is zero in the case of lpTMDPDF, which provides many simplifications. In this notation, the expression for the renormalized lpTMDPDF in eq. (2.17) on a parton reads The expressions for Z TMD g is given in ref. [16], while the expression for the soft factor in δ-regularization is in ref. [25].

JHEP11(2019)121
Given the values of parton matrix elements we find the coefficient functions matching left-and right-hand sides of where f 1,f ←f is the renormalized parton matrix element for PDF operator eq. (2.11), and ⊗ is the short-hand notation for Mellin convolution integral eq. (2.10). Such a relation is valid since the OPE is an operator relation and it is independent of states.
To solve the matching in eq. (3.6) we need the expression for the collinear matrix elements f 1,f ←f . This calculation is trivial in the actual scheme since there is no Lorenzinvariant scale inside the integrands and all loop-integrals for f 1,f ←f are zero in dimensional regularization. For this reason the loop-corrections to f 1,f ←f are given by UV renormalization constant only: where P [1] is the DGLAP evolution kernel at LO. Denoting the perturbative terms for the matching coefficient as we find from eq. (3.6), (3.7) , f ←f (x). (3.10) This procedure cancels the collinear poles that are present in the parton matrix elements. Note that, the last term in eq. (3.10) requires the evaluation of δ L C [1] to order ∼ .

Logarithmic part of the coefficient function
The renormalization group equation allows us to write down the coefficients that accompany the scaling logarithms in the coefficient function. We recall that the coefficient function depends on three scales see eq. (2.10): µ and ζ that are inherited from the TMDPDF, andμ that is the scale of OPE. The behavior on scales µ and ζ is dictated by the TMD evolution equations (2.20), (2.21), while the dependence on scaleμ is canceled by the corresponding dependence of f 1 (x,μ). The latter is given by the DGLAP equation x dy y P f ←f x y , µ f 1,f ←h (y, µ) .

JHEP11(2019)121
Therefore, at the point µ =μ the coefficient function satisfies the pair of equations The solution at NNLO has the simple form The coefficients of logarithms are where Γ g 0 = 4C A is LO cusp anomalous dimension, β 0 = 11/3C A −2/3N f is LO β−function, and we have used that γ g[1] V = −2β 0 . The explicit expressions for these coefficients are given in the appendix B for completeness. The finite parts δ L C (n;0,0) are presented in the next section.
In the expressions above we have setμ = µ, which is a poor choice. In particular, due to this choice one obtains the double-logarithms in the coefficient function and, as the result, a badly convergent perturbative series. A much better behaved coefficient function can be achieved by distinguishing the scales of evolution and OPE. For example, this is realized by applying the ζ-prescription, which consists in the selection of TMD evolution scales along the null-evolution line in the plane (µ, ζ). This line is parameterized as ζ = ζ µ (b), and it is defined by the boundary condition that it passes through the saddle point of the evolution potential [19]. The expression for the coefficient function can be obtained by the substitution (here for gluon distributions) in ζ-prescription: The higher order terms and the derivation of this expression can be found in refs. [19,21]. The coefficient function in ζ-prescription satisfies DGLAP equation, and thus the remaining scale is the OPE scaleμ. In other words, we have

JHEP11(2019)121
where the logarithmic part has simple form (3.19) The finite part δ L C (2,0,0) g←f (x) remains unaffected. Note that, generally the ζ-prescription also modifies the finite part of NNLO expression, as it is happens e.g. for the unpolarized TMDPDF.

Finite part of the coefficient function
In this section we present the finite parts of coefficient function δ L C. The NLO expression read  [9,12,17]. The full -dependent NLO expressions are presented in [17]. The NNLO expressions are where N f is the number of active quark flavors. These expressions is the main result of this work. These results have been recently obtained in ref. [55] with an independent calculation using the exponential regulator of ref. [56] to regularize rapidity divergences. We find full agreement with the final results presented.

lpTMDPDF at NNLO and its contribution Higgs production
The lpTMDPDF and the unpolarized gluon TMDPDF use to be present at the same time in many processes. A particularly important place to study the effect of lpTMDPDF is JHEP11(2019)121 the Higgs production in hadron-hadron collision. In this case the dominating channel for Higgs production is gluon-gluon fusion via the top-quark loop [1], which can be written via an effective interaction term in the Lagrangian [57] where H is the Higgs field, F µν is the gluon field strength tensor, and v is the Higgs vacuum expectation value. The effective coupling constant at NNLO is derived in [58,59].
Using the effective vertex in eq. (4.1) one can derive the TMD factorization theorem for Higgs production following the same steps as in the Drell-Yan case (see refs. [60][61][62][63]). The resulting expression is where y is the Higgs rapidity and x 1,2 = (m 2 H + q 2 T )/se ±y . The function C H is the gluon scalar form-factor (the NNLO expression can be found in [50]), U is the "π 2 -resummation" exponent [64] and the TMD distributions Φ µν are defined in eq. (2.1). For a more accurate and detailed definition we refer to ref. [61]. The scale µ is of the order of the hard scale, m H in this case, and ζ 1 ζ 2 = m 4 H . With the decomposition in eq. (2.3) the product of TMD distributions turns into Therefore, for a consistent phenomenological application of this formula one should consider f 1 and h ⊥ 1 at the same perturbative order. The perturbative inputs up to NNLO are reported in table 1.
It is interesting to mention that if the Higgs boson were a pseudo-scalar particle, then the main change in the structure of cross-section in eq. (4.2) would be a sign of h ⊥ 1 h ⊥ 1 term in eq. (4.3). In this case, the expressions for perturbative corrections in C t and C H are also changed although their LO remains the same [11]. In order to study the numerical impact of our result, the NLO and NNLO matching for lpTMDPDF together with the cross-section in eq. (4.2) have been added to artemide [44,45]. The non-perturbative parts of gluon TMD distributions and gluon rapidity anomalous dimension are unknown, and nowadays the data are not sufficient to fix it. In order to provide some value for a cross section we use the inputs in eq. (2.15)-

JHEP11(2019)121
1N P , where f N P is the non-perturbative function for quarks extracted from a fit of Drell-Yan and Z-boson production data using artemide2.01. The details of this fit have been illustrated in refs. [21,22], and this version of the code takes into account the improvements coming from ref. [66]. The TMD evolution kernel for gluons should be also provided by a non-perturbative part at large value of b, whose precise analytical form is given in [22]. The perturbative calculable parts of the evolution kernel differ in quark and gluon case (at the order that we work) by the Casimir scaling factor C A /C F . Here we have assumed the same scaling for the un-calculable non-perturbative pieces of the evolution kernel. The error band of our prediction come from scale variations of a factor of 2, consistently with ζ-prescription [19].
In order to check the viability of the model assumptions we have compared the cross section in eq. (4.2), integrated in rapidity, with PYTHIA [67,68]. The agreement of our prediction at NNLO and PYTHIA is shown in figure 1 and it is extremely good in the range of q T where the TMD factorization theorem is expected to hold. In that figure we have also included the error provided by PYTHIA, although it is not clearly visible and we have not used any normalization factor.
In figure 2 (left) we have plotted lpTMDPDF, eq. (2.15)-(2.16), as a function of b at x = 0.01 at NLO and at NNLO. The NNLO includes the perturbative correction to the first non-trivial order (which is NLO). This correction appears to be large, say almost a factor 2. The bands show the sensitivity of the distribution to the change of the OPE scalẽ µ → c 4μ with c 4 ∈ (0.5, 2), see eq. (3.18). The relative size of the band decreases between NLO and NNLO. Altogether, this figure points to the fact that the lpTMDPDF effects could have been underestimated up to now. The experimental data on the Higgs differential cross section are still affected by big errors. For a demonstration we have considered the cross section in eq. (4.2) measured at CMS collaboration, where the rapidity is integrated in the interval indicated by that experiment [69]. Because the experimental cross section just uses the data from one particular decay of the Higgs boson we have normalized our cross section with the experimental one integrating in the interval of transverse momenta shown in figure 2 (right). From this figure it is clear that currently the data are not sensitive to the TMD structures.

JHEP11(2019)121
In the Higgs production cross-section the lpTMDPDF mainly affects the low-q T region, as it is demonstrated in figure 3. Practically, the lpTMDPDF can be distinguished from the unpolarized TMDPDF at q T 5-8GeV, where it modifies the values of cross-section by about 5%. Such value of variation band is typical for NNLO approximation, see e.g. [20]. In figure 3 (right) we compare the NNLO cross sections the size of the variation band, which is the maximum deviation value obtained from the variation of all three scales (in ζ-prescription) by factors c i ∈ (0.5, 2) [19]. The variation band is of the order of few percents and the main contribution to it is the µ-band (the scale between hard part and the TMD-evolution factor). Nowadays, these factors can be pushed to N 3 LO reducing the variation band further, if necessary.
Finally, we comment on the positivity relation formulated in ref. [13]:  b-space). Outside of this point the inequality in eq. (4.4) is respected. The situation is exemplified in figure 4, where we plot the ratio of |f 1 |/|h ⊥ 1 | at different values of q T with fixed x (left) and viceversa (right). We also note that the positions of zeros in TMDPDFs strongly depends on the non-perturbative input. In particular, selecting some appropriate model one can, possibly, remove the zero from unpolarized TMDPDF, or fix positions of zeros equal in both gluon TMDPDFs. In other words eq. (4.4) can be used as a serious constraint on non-perturbative part of the TMD distributions. However, we do not see enough theoretical justification for such an approach at the moment.
We have also observed that the ratio |h ⊥ 1 |/|f 1 | tends to saturate at smaller values of x as it is suggested f.i. by [43]. Then for extreme small values of x ∼ 10 −4 it is violated again. However, such values can be outside the applicability region of our calculation since the perturbative expressions for f 1 [16] and h ⊥ 1 (3.22), (3.23) have contributions ∼ a n+1 s ln n (x)/x that should be resummed for a proper comparison.

Conclusions
The gluon transverse momentum dependent parton distribution function (lpTMDPDF) typically accompanies unpolarized gluon TMDPDF within a TMD factorized cross-section. A good example is the factorization formula for the Higgs-production cross-section, where these distributions enter in a plain sum. For this reason, both distributions should be considered at the same order of perturbative accuracy. We have calculated the a 2 s -part (NNLO) for the matching coefficient of lpTMDPDF to twist-2 collinear distributions, which is the main result of this paper. Thanks to this calculation, lpTMDPDF can be considered at the same level of theoretical accuracy as the unpolarized gluon TMDPDF [15,16]. The corresponding formulas are collected in section 3.  for the numerical evaluation of lpTMDPDF is added to the artemide package that can be downloaded from [44,45].
The impact of NNLO correction for lpTMDPDF is very significant and practically doubles the value of the function for moderate b. This fact should not be considered much surprising given that LO term (a 0 s -term) for lpTMDPDF vanishes and the correction that we provide is the one to the first non-null order. The relevance of this effect in the Higgs cross section has been discussed in section 4 and it is resumed in figures 2-3. Unfortunately, at the moment we have not a reliable model for the non-perturbative part of the gluon TMD distribution, and in this work, we have adapted values for distributions extracted in refs. [21,22]. A more detailed study on the non-perturbative part of the gluon TMDPDF is certainly worth in the future. Surprisingly, the model built by us agrees with PYTHIA prediction for low q T values, which are the relevant ones for TMD studies.
In several papers, it has been suggested that unpolarized and linearly polarized gluon TMDPDFs can be measured in association with heavy-quark production [36][37][38][39][40][41][42]. We leave an analysis of these processes for future work because at the moment we miss a full factorization theorem for each of these cases. Nevertheless, the consistency of data with the factorization hypothesis can always be checked with the result provided in this work.

JHEP11(2019)121
A Relevant set of master integrals for linearly polarized gluon TMD Three different types of diagrams arise in the calculation of the unsubstracted TMDPDF matrix element for linearly polarized gluons and the can be addressed on the basis that the exchanged gluons are pure-virtual, virtual-real or real-real. The pure-virtual diagrams, are zero in the dimensional regularization due to the absence of a Lorentz-invariant scale in our scheme of calculation. The virtual-real and real-real diagrams have respectively one and two cut propagators and should be computed directly. The calculation of these two types of diagrams is analogous to the calculation made in refs. [16,24] for the case of unpolarized TMDPDF, albeit with a different Lorentz structure. The main difference and difficulty comes from the term proportional to b µ b ν . The contraction of this term with the projectors generates terms in the numerator as (bq) 2 (where q is a loop-momentum), making the evaluation of the diagrams involved.
For virtual-real diagrams this difficulty can be by-passed by calculating separately virtual subdiagrams. This approach allows to contract the projector only with the real loopmomentum, simplifying the calculation of integrals. For real-real integrals no subdiagrams can be calculated. A set of master integrals in which these diagrams can be decomposed was developed in [16]. In this appendix we present the decomposition of the master integrals original for this work.
A general master integral can be written as where R = {1, (kb) 2 , (kb)(lb), (lb) 2 }. The bold font denotes the scalar product of transverse components only with Euclidian metric. The components k + and l + can be integrated with the help of the introduction of a delta function and they do not enter in the loop-integration (indicated by a d − 1 integral).
The integrals with R = 1, F abcd [1] ≡ F abcd are presented in the appendix C of [16]. In that case, the sum of the indices abcd of the integral is 2. In the present calculation, the new integrals with R = 1 and the sum of the indices abcd is 3. Some of the new integrals can be expressed as a combination of older results, (1−η) 2 (1+x−η) 2 (1−x)F 0020 , (A.5) where B = b 2 /4.
Additionally, we have met three integrals that could not be reduced to a combination of known results: F 1110 [(kb) 2 ], F 1110 [(kb)(lb)], F 1110 [(lb) 2 ]. For these integrals we have derived the expressions in the Schwinger parameterization, and evaluated them in -expansion up to the finite term following the strategy described in the book [70].

JHEP11(2019)121 B Logarithm terms of matching coefficient for lpTMDPDF
In this appendix the logarithmic part of the matching coefficients for lpTMDPDFs are collected. Note that these coefficients are not original, in the sense that they can be predicted from the NLO matching derived in [9,17] via evolution equations as it is described in section 3.2. In our calculation we have derived these expressions directly, as part of the checks.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.