J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Published for SISSA by Springer Received: July 23, 2019 Revised: October 28, 2019 Accepted: November 11, 2019 Published: November 21, 2019 Linearly polarized gluons at next-to-next-to leading order and the Higgs transverse momentum distribution Daniel Gutierrez-Reyes,a Sergio Leal-Gomez,a,b Ignazio Scimemia and Alexey Vladimirovc aDepartamento de F́ısica Teórica and IPARCOS, Universidad Complutense de Madrid (UCM), 28040 Madrid, Spain bFaculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Wien, Austria cInstitut für Theoretische Physik, Universität Regensburg, D-93040 Regensburg, Germany E-mail: dangut01@ucm.es, sergiol95@univie.ac.at, ignazios@ucm.es, alexey.vladimirov@ur.de Abstract: We calculate the small-b (or large-qT ) matching of transverse momentum de- pendent (TMD) distribution for linearly polarized gluons to the integrated gluon distribu- tions 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. Keywords: NLO Computations, QCD Phenomenology ArXiv ePrint: 1907.03780 Open Access, c© The Authors. Article funded by SCOAP3. https://doi.org/10.1007/JHEP11(2019)121 mailto:dangut01@ucm.es mailto:sergiol95@univie.ac.at mailto:ignazios@ucm.es mailto:alexey.vladimirov@ur.de https://arxiv.org/abs/1907.03780 https://doi.org/10.1007/JHEP11(2019)121 J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Contents 1 Introduction 1 2 Gluon TMD distributions 3 2.1 Definition 3 2.2 OPE at small-b 4 2.3 Renormalization of TMDPDF 5 3 Matching coefficient for lpTMDPDF at NNLO 7 3.1 Evaluation of the matching coefficient 7 3.2 Logarithmic part of the coefficient function 9 3.3 Finite part of the coefficient function 11 4 lpTMDPDF at NNLO and its contribution Higgs production 11 5 Conclusions 15 A Relevant set of master integrals for linearly polarized gluon TMD 17 B Logarithm terms of matching coefficient for lpTMDPDF 19 1 Introduction The gluon-gluon fusion is the leading channel for the Higgs boson production in hadron- hadron collisions [1–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–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–12]. Schematically, it reads dσ dyd2qT = σgg→H (2π)2 ∫ d2b e−i(bqT ) ( f1,g(xA, b)f1,g(xB, b) + h⊥1,g(xA, b)h⊥1,g(xB, b) ) , (1.1) where σgg→H is the factorized gluon-gluon-Higgs cross-section, xA,B are the collinear frac- tions of gluon momenta, f1 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]. – 1 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 In TMD factorization each TMD distribution (f1,g and h⊥1,g in this case) is an inde- pendent 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 coeffi- cient 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–18]. Such a high order is required because of the sizes of theo- retical 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 high- and low-energy data [21, 22]. The small-b limit of the unpolarized gluon TMDPDF, f1, 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 esti- mating 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–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 ap- proach. In fact, using standard resummation, see e.g. [12, 20, 27–30], the small-b expan- sion 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 qT , however one has different requirements in the realization of the perturbative series. So, while in the usual resumma- tion 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–35]. Among these it is worth a special mentioning the case of heavy-quark production [36–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], |h⊥1,g←h(x, qT )|/|f1,g←h(x, qT )| ≤ 1. (1.2) 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. 1In ref. [23] the authors use the differential cross section for Higgs production at NNLO which includes only the NLO matching coefficient for the linearly polarized gluons. This counting is different from the one required by TMD factorization as explained in the text. – 2 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 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-N3LL. The results of the phenomenological analysis are discussed in section 4. 2 Gluon TMD distributions 2.1 Definition The TMD distribution of gluons in a hadron is given by the following matrix element Φg←h,µν(x, b) = 1 xp+ ∫ dλ 2π e−ixp +λ (2.1) × 〈P, S|T̄ { F+µ (λn+ b) W̃n (λn+ b) } T { W̃ †n(0)F+ν(0) } |P, S〉, where n is a light-like vector, Fµν is the gluon field strength tensor, and W̃ denotes the half-infinite Wilson line in the direction n W̃n(z) = P exp ( ig ∫ 0 −∞ dσA+(nσ + z) ) . (2.2) The Wilson lines W̃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 n2 = 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 Φµν g←h(x, b) = − gµνT 2(1− ε) f1,g←h(x, b) + h⊥1,g←h(x, b) ( gµνT 2(1− ε) + bµbν b2 ) , (2.3) where b2 = −b2 > 0. For future necessity, the decomposition in eq. (2.3) is given in d = 4−2ε-dimensions as it was defined in [15, 17]. Both f1 and h⊥1 contribute to the gluon- induced 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] – 3 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 (here dots denote the staple gauge link, as in (2.1)), f1,g←h(x, b) = −gµνT Φg←h,µν(x, b) = 1 xp+ ∫ dλ 2π e−ixp +λ〈P |F+µ (λn+ b) . . . F+µ(0)|P 〉, (2.4) whereas the linearly-polarized tensor is orthogonal to it. In turn the lpTMDPDF is given by h⊥1,g←h(x, b) = 1 1− 2ε ( gµνT + 2(1− ε)b µbν b2 ) Φg←h,µν(x, b). (2.5) 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), Φg←h,µν(x,k) = ∫ d2b (2π)2 ei(bk)Φg←h,µν(x, b) (2.6) = − gµνT 2 f1,g←h(x,k) + h⊥1,g←h(x,k) ( gµνT 2 + kµkν k2 ) , where f1,g←h(x,k) = ∫ ∞ 0 |b|d|b| 2π J0(|b||k|) f1,g←h(x, b), (2.7) h⊥1,g←h(x,k) = − ∫ ∞ 0 |b|d|b| 2π J2(|b||k|)h⊥1,g←h(x, b). (2.8) 2.2 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 f1,g←h(x, b;µ, ζ) = ∑ f ∫ 1 x dy y Cg←f (y, b;µ, ζ; µ̃) f1,f←h ( x y , µ̃ ) +O(b2) (2.9) h⊥1,g←h(x, b;µ, ζ) = ∑ f ∫ 1 x dy y δLCg←f (y, b;µ, ζ; µ̃) f1,f←h ( x y , µ̃ ) +O(b2), (2.10) where the sum runs over the active parton flavors (quarks and gluon), and f1(x, µ) is unpolarized collinear distributions defined as usual f1,q←h(x, µ) = ∫ dλ 2π e−ixp +λ〈P |T̄{q̄ (λn) W̃n(λn)}γ + 2 T{W̃ †n(0)q(0)}|P 〉, (2.11) f1,g←h(x, µ) = 1 xp+ ∫ dλ 2π e−ixp +λ〈P |T̄ { F+µ(λn)W̃n(λn) } T { W̃ †n(0)F+µ(0) } |P 〉. (2.12) – 4 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Concerning the notation, here and in the following we distinguish the unpolarized TMD- PDF f1(x, b) and unpolarized collinear PDF f1(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 δLC, are to be calculated in QCD perturbation theory. The three-order calculation yields Cg←f (x, b;µ, ζ; µ̃) = δgfδ(1− x) +O(as), (2.13) δLCg←f (x, b;µ, ζ; µ̃) = O(as), (2.14) where as = g2/(4π)2 is QCD coupling constant. Nowadays, the coefficients Cf←h(x, b) are known at a2 s-order (NNLO) [15, 16, 24, 47], whereas coefficients δLCf←h(x, b) are known at as-order (NLO2) [9, 12, 17]. In the following section we present NNLO expression for δLCg←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 b2 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 h⊥1,g←h(x, b) = ∑ f ∫ 1 x dy y δLCg←f (y, b) f1,f←h ( x y ) h⊥1NP(x, y, b2), (2.15) and a similar expression can be used for f1(x, b) with a different f1NP(x, y, b2) 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 lim b2→0 h⊥1NP(x, y, b2) ' 1 +O(b2), (2.16) 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 fNP 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 fNP = f1NP = h⊥1NP . We comment about the consistency of this choice in section 4. 2.3 Renormalization of TMDPDF The TMD operator contains ultraviolet (UV) and rapidity divergences. Both these diver- gences 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) 2In literature related to TMD calculations, e.g. in refs. [9, 17], the orders of δLCf←h are traditionally counted alike the unpolarized case. So, the linear as-terms are denoted as NLO. Here we use the same convention. – 5 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 and ζ (the rapidity divergences renormalization scale). The renormalized expression for the TMD distribution Φg←h reads Φg←h(x, b;µ, ζ) = ZTMD g (µ, ζ|ε)Rg ( b, µ, ζ|ε, δ + p+ ) Φunsub. g←h ( x, b|ε, δ + p+ ) , (2.17) where Φµν;unsub. g←h denotes the bare or unsubtracted TMD distribution, either f1 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 di- mensional regularization (d = 4−2ε) that regularizes UV divergences, that are renormalized by the factor Zg; δ is the parameter of δ-regularization [16, 25] which regularizes rapidity divergences that are renormalized by the factor Rg. The renormalization factors Zg and Rg 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], S(b) = Trcolor Nc 〈0| [ W T † n W̃ T n̄ ] (b) [ W̃ T † n̄ W T n ] (0)|0〉, (2.18) where W̃n and W̃n̄ are Wilson lines along n and n̄ (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] Rg ( b, µ, ζ|ε, δ + p+ ) = S−1/2 ( b|ε, δ = δ+ 2p+ √ ζ ) . (2.19) The variable p+ is parton momentum [46], and is required to define the Lorentz invari- ant 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 con- dition 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 µ2 d dµ2 Φg←h(x, b;µ, ζ) = γg(µ, ζ) 2 Φg←h(x, b;µ, ζ), (2.20) ζ d dζ Φg←h(x, b;µ, ζ) = −Dg(µ, b)Φg←h(x, b;µ, ζ). (2.21) These equations are the same for all gluon TMD distributions of leading twist. The anoma- lous dimensions are defined via the corresponding renormalization constants and they are – 6 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 known up to three-loop order inclusively [49–52]. Note, that the renormalization factor Zg also contains the gluon-field renormalization part, therefore, γG = 2ÂD(Z3 − Zg) (2.22) where the symbol ÂD extracts the coefficient of ε−1 with a pre-factor n! at the nth pertur- bative order. Anomalous dimensions γg and D satisfy the integrability condition (also known as Collins-Soper equation [53]) 2µ2dDg(b, µ) dµ2 = −ζ dγ g(µ, ζ) dζ = Γgcusp(µ), (2.23) 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 γg(µ, ζ) = Γgcusp(µ) ln ( µ2 ζ ) − γgV , (2.24) where γgV is anomalous dimension of the vector form factor for gluon. The rapidity anoma- lous dimension Dg has not such a simple representation due to the presence of an extra dimensional parameter b2. It generally contains all powers of logarithms ln(µ2b2), that at some large values of b2 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: Φg←h(x, b;µ1, ζ1) = Rg[b; (µ1, ζ1)→ (µ2, ζ2)]Φg←h(x, b;µ2, ζ2), (2.25) where the TMD renormalization factor reads Rg[b; (µ1, ζ1)→ (µ2, ζ2)] = exp [∫ P ( γg(µ, ζ) dµ µ −Dg(µ, b) dζ ζ )] . (2.26) 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 f1 and lpTMDPDF h⊥1 . 3 Matching coefficient for lpTMDPDF at NNLO 3.1 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: – 7 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 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 diag. = (b2)2ε ( g1(x, ε) + ( δ+ p+ )ε g2(x, ε) + ( δ+ p+ )−ε g3(x, ε) (3.1) + ln ( δ+ p+ ) g4(x, ε) + ln2 ( δ+ p+ ) g5(x, ε) ) . The functions g2 and g3 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 g4 and g5 are canceled by the rapidity renormalization factor. However, due to the absence of three-order term, the functions g5 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 h⊥;unsub. 1;f←f ′ (x, b) = Φunsub. f←f ′ (x, b) = ∞∑ n=1 ansΦ [n]unsub. f←f ′ , S(b) = 1 + ∞∑ n=1 ansS [n], (3.2) where as = g2/(4π)2. The tree-order term is zero in the case of lpTMDPDF, Φ [0]unsub. f←f ′ = 0, (3.3) which provides many simplifications. In this notation, the expression for the renormalized lpTMDPDF in eq. (2.17) on a parton reads Φ [1] f←f ′ = Φ [1]unsub. f←f ′ (3.4) Φ [2] f←f ′ = Φ [2]unsub. f←f ′ − S[1]Φ [1]unsub. f←f ′ 2 + Z [1]TMD g Φ [1]unsub. f←f ′ . (3.5) The expressions for ZTMD g is given in ref. [16], while the expression for the soft factor in δ-regularization is in ref. [25]. – 8 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Given the values of parton matrix elements we find the coefficient functions matching left- and right-hand sides of h⊥1,g←f (x, b) = ∑ f=g,q,q̄ [ δLCg←f ′(b)⊗ f1,f ′←f ] (x), (3.6) where f1,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 f1,f ′←f . This calculation is trivial in the actual scheme since there is no Lorenz- invariant scale inside the integrands and all loop-integrals for f1,f ′←f are zero in dimensional regularization. For this reason the loop-corrections to f1,f ′←f are given by UV renormal- ization constant only: f [0] 1,f←f ′(x) = δff ′δ(1− x), f [1] 1,f←f ′(x) = − P [1] f←f ′(x) ε , (3.7) where P [1] is the DGLAP evolution kernel at LO. Denoting the perturbative terms for the matching coefficient as δLCg←f (x, b) = ∞∑ n=1 ans δ LC [n] g←f (x, b), (3.8) we find from eq. (3.6), (3.7), δLC [0] g←f (x, b) = 0, δLC [1] g←f (x, b) = h ⊥[1] g←f ′(x, b), (3.9) δLC [2] g←f (x, b) = h ⊥[2] g←f ′(x, b) + 1 ε ∑ f ′ [ δLC [1] g←f ′(b)⊗ P [1] 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 δLC [1] to order ∼ ε. 3.2 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 f1(x, µ̃). The latter is given by the DGLAP equation µ2 d dµ2 f1,f←h(x, µ) = ∑ f ′=g,q,q̄ ∫ 1 x dy y Pf←f ′ ( x y , µ ) f1,f ′←h(y, µ) . (3.11) – 9 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Therefore, at the point µ = µ̃ the coefficient function satisfies the pair of equations µ2 d dµ2 δLCg←f (x, b;µ, ζ, µ) (3.12) = ∑ f ′=g,q,q̄ ∫ 1 x dy y δLCg←f ′ ( x y , b;µ, ζ, µ )( γgV (µ, ζ) 2 δff ′δ(ȳ)− Pf ′←f (y, µ) ) , ζ d dζ δLCg←f (x, b;µ, ζ, µ) = −Dg(µ, b)δLCg←f (x, b;µ, ζ, µ). (3.13) The solution at NNLO has the simple form δLC [2] g←f (x, b;µ, ζ, µ) = ( −1 2 L2 µ + Lµlζ ) δLC (2,1,1) g←f (x) + Lµδ LC (2,1,0) g←f (x) + δLC (2,0,0) g←f (x), (3.14) where Lµ = ln ( µ2b2 4e−γE ) , lζ = ln ( µ2 ζ ) . (3.15) The coefficients of logarithms are δLC (2,1,1) g←f (x) = Γg0 2 δLC (1,0,0) g←f (x), (3.16) δLC (2,1,0) g←f (x) = 2β0δ LC (1,0,0) g←f (x)− ∑ f ′=g,q,q̄ [δLC (1,0,0) g←f ′ ⊗ P [1] f ′←f ](x), where Γg0 = 4CA is LO cusp anomalous dimension, β0 = 11/3CA−2/3Nf 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 δLC(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: lζ = Lµ 2 − 2β0 Γg0 +O(as). (3.17) 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 δLCg←f (x, b;µ, ζµ(b), µ̃) = δLCg←f (x, b; µ̃), (3.18) – 10 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 where the logarithmic part has simple form δLC [2] g←f (x, b; µ̃) = β0δ LC (1,0,0) g←f (x)− ∑ f ′=g,q,q̄ [ δLC (1,0,0) g←f ′ ⊗ P [1] f ′←f ] (x) Lµ̃ + δLC (2,0,0) g←f (x). (3.19) The finite part δLC (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. 3.3 Finite part of the coefficient function In this section we present the finite parts of coefficient function δLC. The NLO expression read δLC(1,0,0) g←g (x, b) = −CA 4(1− x) x , (3.20) δLC(1,0,0) g←q (x, b) = −CF 4(1− x) x , (3.21) where CA = Nc(= 3) and CF = (N2 c − 1)/2Nc(= 4/3) are eigenvalues of quadratic Casimir operators for adjoint and fundamental representations in SU(Nc)(SU(3)) group. The result in eq. (3.20), (3.21) agrees with [9, 12, 17]. The full ε-dependent NLO expressions are presented in [17]. The NNLO expressions are δLC(2;0,0) g←g (x) =C2 A [ −16 1−x x (Li2(x)−lnx)+ 124 3 lnx+ ( 148 9 +20ζ2 ) 1−x x −8ln2x− 100 9 (1−x)− 4 9 x(11x−14) ] +CFNf ·4 [ ln2x−2 (1−x)3 x ] +CANf · 4 9 [ 17 1−x x +1−3x−x2+6lnx ] , (3.22) δLC(2;0,0) g←q (x) =CF (CF−CA) [ 8 1−x x (ln(1−x)+ln2(1−x))−20lnx+4ln2x+8(1−x) ] +CFCA [ 16 1−x x ( 11 18 + 5 4 ζ2− ln(1−x) 3 −Li2(x) ) +4 lnx x (4+5x−x lnx) ] +CFNf · 16 9 1−x x [2+3ln(1−x)], (3.23) where Nf 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. 4 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 – 11 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Function H Cg←f , δLCg←f Γcusp D γF αs running PDF evolution NLO αs αs α2 s αs resummed α2 s NLO provided by NNPDF3.1 [65] NNLO α2 s α2 s α3 s α2 s resummed α3 s NNLO provided by NNPDF3.1 [65] Table 1. Summary of perturbative orders used for each part of the cross section. The symbol H stands for the first line of eq. (4.2). 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] LggH = as(µ)Ct(µ) 3v FAµνF A,µνH, (4.1) 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–63]). The resulting expression is dσ dyd2qT = 2σ0(µ) π C2 t (µ)U(µ,−µ)|CH(−m2 H ,−µ2)|2 (4.2)∫ d2b 4π ei(bqT )Φµν g←h1(x1, b;µ, ζ1)Φµν g←h2(x2, b;µ, ζ2), where y is the Higgs rapidity and x1,2 = √ (m2 H + q2 T )/se±y. The function CH 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, mH in this case, and ζ1ζ2 = m4 H . With the decomposition in eq. (2.3) the product of TMD distributions turns into Φµν g←h1(x1,b)Φµν g←h2(x2,b) = 1 2 ( f1,g←h1(x1,b)f1,g←h2(x2,b)+h⊥1,g←h1(x1,b)h⊥1,g←h2(x2,b) ) . (4.3) Therefore, for a consistent phenomenological application of this formula one should consider f1 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 Ct and CH are also changed although their LO remains the same [11]. – 12 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Figure 1. The cross section in eq. (4.2) integrated over all rapidity range with artemide2.01 at NNLO and PYTHIA. The errors of PYTHIA are included, although not clearly visible. The shaded area shows the variation band in µ̃, see eq. (3.18). In order to study the numerical impact of our result, the NLO and NNLO match- ing 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)– (2.16) with fNP = f1NP = h⊥1NP , where fNP 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 glu- ons 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 CA/CF . 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 qT 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 scale µ̃→ c4µ̃ with c4 ∈ (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. – 13 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Figure 2. (left) The lpTMDPDF, eq. (2.15)–(2.16), as a function of b at x = 0.01. The shaded area shows the variation band in µ̃, see eq. (3.18). (right) Comparison of Higgs-production cross-section with variation band to the measurement presented in [69] by CMS collaboration. 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 partic- ular 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. In the Higgs production cross-section the lpTMDPDF mainly affects the low-qT region, as it is demonstrated in figure 3. Practically, the lpTMDPDF can be distinguished from the unpolarized TMDPDF at qT . 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 ci ∈ (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 N3LO reducing the variation band further, if necessary. Finally, we comment on the positivity relation formulated in ref. [13]: |f1(x, qT )| − |h⊥1 (x, qT )| > 0. (4.4) This relation is a consequence of positive definiteness of the gluon-polarization matrix in a free theory, and certainty hold at LO. However, it does not need to be accomplished at higher order in perturbation theory. The positivity bound is formulated in momentum space, whereas all perturbative calculation are performed in coordinate space. This causes an additional problem since the Hankel transform of a positive function is not necessary a positive function. Within our model we have checked that it is easy to get a violation of this bound, for any fixed value of x and qT . Typically, the violation happens in the vicinity of sign change point of f1 (note, that our realization of f1 is positive-definite in – 14 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Figure 3. Cross section for Higgs production including linearly polarized gluon effects at different orders. (left) The motion of center lines of cross-section integrated over all rapidities at different perturbative orders for lpTMDPDF. The black (blue) lines correspond the case of positive (nega- tive) contribution for h⊥1 h ⊥ 1 -term in eq. (4.3). (right) The scale-variation band for the cross section at NNLO. The parity even (odd) Higgs case is represented with a green (orange) band. In both figures the center of mass energy is set as in [69] and rapidity is integrated over its complete range. 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 |f1|/|h⊥1 | at different values of qT 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 |/|f1| 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 f1 [16] and h⊥1 (3.22), (3.23) have contributions ∼ an+1 s lnn(x)/x that should be resummed for a proper comparison. 5 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 a2 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.2, 3.3. They are also attached to the publication as supplementary material in the form of Mathematica-notebook. The module – 15 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 Figure 4. Ratio of linearly polarized and unpolarized gTMD to check eq. (1.2) as a function of qT at fixed x = 0.01 (left) and as a function of x at fixed qT = 1 GeV (right). 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 (a0 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 qT 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–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. Acknowledgments D.G.R., S.L.G. and I.S. are supported by the Spanish MECD grant FPA2016-75654-C2-2- P. This project has received funding from the European Union Horizon 2020 research and innovation program under grant agreement No 824093 (STRONG-2020). D.G.R. acknowl- edges the support of the Universidad Complutense de Madrid through the predoctoral grant CT17/17-CT18/17. S.L.G. is supported by the Austrian Science Fund FWF under the Doctoral Program W1252-N27 Particles and Interactions. – 16 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 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 loop- momentum, 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 Fabcd[R] = (2π)2 ∫ dd−1k dd−1l (2π)2d Rei(kb)ei(lb)δ(k2)θ(−k−)δ(l2)θ(−l−) [(l + p)2]a[(k + p)2]b[(k + l + p)2]c[(k + l)2]d , (A.1) where R = {1, (kb)2, (kb)(lb), (lb)2}. The bold font denotes the scalar product of trans- verse components only with Euclidian metric. The components k+ and l+ can be integrated with the help of the introduction of a delta function 1 = ∫ ∞ −∞ dη p+δ ( (1− η)p+ + l+ ) (A.2) and they do not enter in the loop-integration (indicated by a d− 1 integral). The integrals with R = 1, Fabcd[1] ≡ Fabcd 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 6= 1 and the sum of the indices abcd is 3. Some of the new integrals can be expressed as a combination of older results, F0210[(kb)2]/B = 2 ( (1+2ε)(x−η)− ε(1−2ε) 1+ε 1−η x ) F0110− 2(1−2ε) 1+ε (1−x)F0020, (A.3) F0210[(kb)(lb)]/B = 2(1−2ε) 1+ε 1−η 1+x−η ( ε 1−η x F0110−(1+ε)(η−x)F0110+(1−x)F0020 ) +2(1+x−η)F(−1)210+2ηF0110, (A.4) – 17 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 F0210[(lb)2]/B = x (1+x−η)2 ( 2(1−2ε) 1+ε (x(η−x)−(1+ε)(1−η))F0110 +ε 2(1−2ε) 1+ε x(1−x)F0020−2(1−2ε)(1−η)2F0110 ) −4(1−η)F(−1)210 − 1 (1+x−η)2 ( 2ε(1−2ε) 1+ε (1−η)3 x −2η(1−2ε)(1−η)2 ) F0110 − 2(1−2ε) 1+ε (1−η)2 (1+x−η)2 (1−x)F0020, (A.5) F0120[(kb)2]/B = ( 4(x−η)+ 2(1−2ε) 1+ε 1−x x (1+x−η) ) F0020 +ε 2(1−2ε) 1+ε 1−η x2 (1+x−η)F0110, (A.6) F0120[(kb)(lb)]/B =−2(1−2ε) 1+ε 1−η x ( ε(1−η) x F0110+(1−x)F0020 ) −2F0110+2(1+x−η)F(−1)210+2ηF0020, (A.7) F0120[(lb)2]/B = 2(1−2ε) 1+ε (1−η)2 x(1+x−η) ( ε 1−η x F0110+(1+ε) x 1−η F0110+(1−x)F0020 ) −4(1−η)F(−1)210, (A.8) F1020[(kb)2]/B = 2(1−2ε) 1+ε (η−x)2 xη ( ε η−x x F1010+(1+ε) x η−x F1010+(1−x)F0020 ) +4(x−η)F1(−1)20, (A.9) F1020[(kb)(lb)]/B =−2(1−2ε) 1+ε η−x x ( ε η−x x F1010+(1−x)F0020 ) −2F1010+2(1+x−η)F0020+2ηF1(−1)20, (A.10) F1020[(lb)2]/B = 2(1−2ε) 1+ε η x ( ε η−x x F1010+(1−x)F0020 ) −4(1−η)F0020, (A.11) F0021[(kb)2]/B =−2(1−2ε) 1−ε η−x 1−x ((1+x−η)F0020−(η−x)F0011) +4(x−η) ( F0011−F0020−F(−1)021 ) , (A.12) F0021[(kb)(lb)]/B = 2(1−2ε) 1−ε (η−x)(1−η) 1−x (F0020+F0011)+2(1+x−2η)F(−1)021 −2(1−η)(F0011−F0020)−2F0020, (A.13) F0021[(lb)2]/B =−2(1−2ε) 1−ε 1−η 1−x (ηF0020−(1−η)F0011)−4(1−η)F(−1)021, (A.14) where B = b2/4. Additionally, we have met three integrals that could not be reduced to a combina- tion of known results: F1110[(kb)2], F1110[(kb)(lb)], F1110[(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]. – 18 – J H E P 1 1 ( 2 0 1 9 ) 1 2 1 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 pre- dicted 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. Recalling that the perturbative expansion of the coefficient function in eq. (2.10) is δLCg←f (x, b;µ, ζ, µ) = ∞∑ n=1 ans δ LC [n] g←f (x, b;µ, ζ, µ), (B.1) with as = g2/(4π)2 and solving the system of eq. (3.12), (3.13) we obtain δLC [1] g←f (x,b;µ,ζ,µ) = δLC (1,0,0) g←f (x), (B.2) δLC [2] g←f (x,b;µ,ζ,µ) = ( −1 2 L2 µ+Lµlζ ) δLC (2,1,1) g←f (x)+Lµδ LC (2,1,0) g←f (x)+δLC (2,0,0) g←f (x), (B.3) where Lµ = ln ( µ2b2 4e−γE ) , lζ = ln ( µ2 ζ ) . (B.4) Using expression for the NLO coefficients (3.20), (3.21) and the LO DGLAP kernels [71] and expressions for anomalous dimensions (see e.g. [16]) we obtain δLC(2,1,1) g←g (x) = −8C2 A 1− x x , (B.5) δLC(2;1,0) g←g (x) = −16C2 A { 1 + x x lnx+ 1− x x [ x 6 (2− x) + 15 4 − ln(1− x) ]} (B.6) + 16CFTrNf [ 1 3 1− x x (2 + (2− x)x) + lnx ] + 16 3 CATrNf 1− x x , δLC(2,1,1) g←q (x) = −8CFCA 1− x x , (B.7) δLC(2;1,0) g←q (x) = −4CFCA [ 1− x x ( 43 3 + x ) + 4 1 + x x lnx ] (B.8) + 4C2 F [ 1− x x (x+ 4 ln(1− x)) + 2 lnx ] + 32 3 CFTrNf 1− x x , where CA = Nc, CF = (N2 c − 1)/2Nc are Casimir eigenvalues of adjoint and fundamen- tal representation for SU(Nc)-gauge group, Tr = 1/2 is the normalization of Gell-Mann matrices, and Nf is the number of quark flavors. 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. – 19 – https://creativecommons.org/licenses/by/4.0/ J H E P 1 1 ( 2 0 1 9 ) 1 2 1 References [1] J.R. Ellis, M.K. Gaillard and D.V. Nanopoulos, A Phenomenological Profile of the Higgs Boson, Nucl. Phys. B 106 (1976) 292 [INSPIRE]. [2] M. Spira, A. Djouadi, D. Graudenz and P.M. Zerwas, Higgs boson production at the LHC, Nucl. Phys. B 453 (1995) 17 [hep-ph/9504378] [INSPIRE]. [3] A. Djouadi, The Anatomy of electro-weak symmetry breaking. I: The Higgs boson in the standard model, Phys. Rept. 457 (2008) 1 [hep-ph/0503172] [INSPIRE]. [4] S. Catani, E. D’Emilio and L. Trentadue, The Gluon Form-factor to Higher Orders: Gluon Gluon Annihilation at Small Q−transverse, Phys. Lett. B 211 (1988) 335 [INSPIRE]. [5] J. Collins, Foundations of perturbative QCD, Cambridge University Press (2013) [INSPIRE]. [6] M.G. Echevarria, A. Idilbi and I. Scimemi, Factorization Theorem For Drell-Yan At Low qT And Transverse Momentum Distributions On-The-Light-Cone, JHEP 07 (2012) 002 [arXiv:1111.4996] [INSPIRE]. [7] M.G. Echevarria, A. Idilbi and I. Scimemi, Unified treatment of the QCD evolution of all (un-)polarized transverse momentum dependent functions: Collins function as a study case, Phys. Rev. D 90 (2014) 014003 [arXiv:1402.0869] [INSPIRE]. [8] A. Vladimirov, Structure of rapidity divergences in multi-parton scattering soft factors, JHEP 04 (2018) 045 [arXiv:1707.07606] [INSPIRE]. [9] M.G. Echevarria, T. Kasemets, P.J. Mulders and C. Pisano, QCD evolution of (un)polarized gluon TMDPDFs and the Higgs qT -distribution, JHEP 07 (2015) 158 [Erratum ibid. 05 (2017) 073] [arXiv:1502.05354] [INSPIRE]. [10] S. Catani and M. Grazzini, QCD transverse-momentum resummation in gluon fusion processes, Nucl. Phys. B 845 (2011) 297 [arXiv:1011.3918] [INSPIRE]. [11] D. Boer, W.J. den Dunnen, C. Pisano, M. Schlegel and W. Vogelsang, Linearly Polarized Gluons and the Higgs Transverse Momentum Distribution, Phys. Rev. Lett. 108 (2012) 032002 [arXiv:1109.1444] [INSPIRE]. [12] T. Becher, M. Neubert and D. Wilhelm, Higgs-Boson Production at Small Transverse Momentum, JHEP 05 (2013) 110 [arXiv:1212.2621] [INSPIRE]. [13] P.J. Mulders and J. Rodrigues, Transverse momentum dependence in gluon distribution and fragmentation functions, Phys. Rev. D 63 (2001) 094021 [hep-ph/0009343] [INSPIRE]. [14] T. Becher and M. Neubert, Drell-Yan Production at Small qT , Transverse Parton Distributions and the Collinear Anomaly, Eur. Phys. J. C 71 (2011) 1665 [arXiv:1007.4005] [INSPIRE]. [15] T. Gehrmann, T. Luebbert and L.L. Yang, Calculation of the transverse parton distribution functions at next-to-next-to-leading order, JHEP 06 (2014) 155 [arXiv:1403.6451] [INSPIRE]. [16] M.G. Echevarria, I. Scimemi and A. Vladimirov, Unpolarized Transverse Momentum Dependent Parton Distribution and Fragmentation Functions at next-to-next-to-leading order, JHEP 09 (2016) 004 [arXiv:1604.07869] [INSPIRE]. [17] D. Gutiérrez-Reyes, I. Scimemi and A.A. Vladimirov, Twist-2 matching of transverse momentum dependent distributions, Phys. Lett. B 769 (2017) 84 [arXiv:1702.06558] [INSPIRE]. – 20 – https://doi.org/10.1016/0550-3213(76)90382-5 https://inspirehep.net/search?p=find+J+%22Nucl.Phys.,B106,292%22 https://doi.org/10.1016/0550-3213(95)00379-7 https://arxiv.org/abs/hep-ph/9504378 https://inspirehep.net/search?p=find+EPRINT+hep-ph/9504378 https://doi.org/10.1016/j.physrep.2007.10.004 https://arxiv.org/abs/hep-ph/0503172 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0503172 https://doi.org/10.1016/0370-2693(88)90912-4 https://inspirehep.net/search?p=find+J+%22Phys.Lett.,B211,335%22 https://inspirehep.net/search?p=find+IRN+9162291 https://doi.org/10.1007/JHEP07(2012)002 https://arxiv.org/abs/1111.4996 https://inspirehep.net/search?p=find+EPRINT+arXiv:1111.4996 https://doi.org/10.1103/PhysRevD.90.014003 https://arxiv.org/abs/1402.0869 https://inspirehep.net/search?p=find+EPRINT+arXiv:1402.0869 https://doi.org/10.1007/JHEP04(2018)045 https://doi.org/10.1007/JHEP04(2018)045 https://arxiv.org/abs/1707.07606 https://inspirehep.net/search?p=find+EPRINT+arXiv:1707.07606 https://doi.org/10.1007/JHEP07(2015)158 https://arxiv.org/abs/1502.05354 https://inspirehep.net/search?p=find+EPRINT+arXiv:1502.05354 https://doi.org/10.1016/j.nuclphysb.2010.12.007 https://arxiv.org/abs/1011.3918 https://inspirehep.net/search?p=find+EPRINT+arXiv:1011.3918 https://doi.org/10.1103/PhysRevLett.108.032002 https://doi.org/10.1103/PhysRevLett.108.032002 https://arxiv.org/abs/1109.1444 https://inspirehep.net/search?p=find+EPRINT+arXiv:1109.1444 https://doi.org/10.1007/JHEP05(2013)110 https://arxiv.org/abs/1212.2621 https://inspirehep.net/search?p=find+EPRINT+arXiv:1212.2621 https://doi.org/10.1103/PhysRevD.63.094021 https://arxiv.org/abs/hep-ph/0009343 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0009343 https://doi.org/10.1140/epjc/s10052-011-1665-7 https://arxiv.org/abs/1007.4005 https://inspirehep.net/search?p=find+EPRINT+arXiv:1007.4005 https://doi.org/10.1007/JHEP06(2014)155 https://arxiv.org/abs/1403.6451 https://inspirehep.net/search?p=find+EPRINT+arXiv:1403.6451 https://doi.org/10.1007/JHEP09(2016)004 https://arxiv.org/abs/1604.07869 https://inspirehep.net/search?p=find+EPRINT+arXiv:1604.07869 https://doi.org/10.1016/j.physletb.2017.03.031 https://arxiv.org/abs/1702.06558 https://inspirehep.net/search?p=find+EPRINT+arXiv:1702.06558 J H E P 1 1 ( 2 0 1 9 ) 1 2 1 [18] M.-X. Luo, X. Wang, X. Xu, L.L. Yang, T.-Z. Yang and H.X. Zhu, Transverse Parton Distribution and Fragmentation Functions at NNLO: the Quark Case, JHEP 10 (2019) 083 [arXiv:1908.03831] [INSPIRE]. [19] I. Scimemi and A. Vladimirov, Systematic analysis of double-scale evolution, JHEP 08 (2018) 003 [arXiv:1803.11089] [INSPIRE]. [20] J. Cruz-Martinez, T. Gehrmann, E.W.N. Glover and A. Huss, Second-order QCD effects in Higgs boson production through vector boson fusion, Phys. Lett. B 781 (2018) 672 [arXiv:1802.02445] [INSPIRE]. [21] I. Scimemi and A. Vladimirov, Analysis of vector boson production within TMD factorization, Eur. Phys. J. C 78 (2018) 89 [arXiv:1706.01473] [INSPIRE]. [22] V. Bertone, I. Scimemi and A. Vladimirov, Extraction of unpolarized quark transverse momentum dependent parton distributions from Drell-Yan/Z-boson production, JHEP 06 (2019) 028 [arXiv:1902.08474] [INSPIRE]. [23] X. Chen et al., Precise QCD Description of the Higgs Boson Transverse Momentum Spectrum, Phys. Lett. B 788 (2019) 425 [arXiv:1805.00736] [INSPIRE]. [24] M.G. Echevarria, I. Scimemi and A. Vladimirov, Transverse momentum dependent fragmentation function at next-to-next-to-leading order, Phys. Rev. D 93 (2016) 011502 [Erratum ibid. D 94 (2016) 099904] [arXiv:1509.06392] [INSPIRE]. [25] M.G. Echevarria, I. Scimemi and A. Vladimirov, Universal transverse momentum dependent soft function at NNLO, Phys. Rev. D 93 (2016) 054004 [arXiv:1511.05590] [INSPIRE]. [26] D. Gutierrez-Reyes, I. Scimemi and A. Vladimirov, Transverse momentum dependent transversely polarized distributions at next-to-next-to-leading-order, JHEP 07 (2018) 172 [arXiv:1805.07243] [INSPIRE]. [27] G. Bozzi, S. Catani, D. de Florian and M. Grazzini, Transverse-momentum resummation and the spectrum of the Higgs boson at the LHC, Nucl. Phys. B 737 (2006) 73 [hep-ph/0508068] [INSPIRE]. [28] S. Mantry and F. Petriello, Factorization and Resummation of Higgs Boson Differential Distributions in Soft-Collinear Effective Theory, Phys. Rev. D 81 (2010) 093007 [arXiv:0911.4135] [INSPIRE]. [29] D. de Florian, G. Ferrera, M. Grazzini and D. Tommasini, Transverse-momentum resummation: Higgs boson production at the Tevatron and the LHC, JHEP 11 (2011) 064 [arXiv:1109.2109] [INSPIRE]. [30] W. Bizoń et al., Fiducial distributions in Higgs and Drell-Yan production at N3LL+NNLO, JHEP 12 (2018) 132 [arXiv:1805.05916] [INSPIRE]. [31] D. Boer, S.J. Brodsky, P.J. Mulders and C. Pisano, Direct Probes of Linearly Polarized Gluons inside Unpolarized Hadrons, Phys. Rev. Lett. 106 (2011) 132001 [arXiv:1011.4225] [INSPIRE]. [32] A. Metz and J. Zhou, Distribution of linearly polarized gluons inside a large nucleus, Phys. Rev. D 84 (2011) 051503 [arXiv:1105.1991] [INSPIRE]. [33] F. Dominguez, J.-W. Qiu, B.-W. Xiao and F. Yuan, On the linearly polarized gluon distributions in the color dipole model, Phys. Rev. D 85 (2012) 045003 [arXiv:1109.6293] [INSPIRE]. – 21 – https://doi.org/10.1007/JHEP10(2019)083 https://arxiv.org/abs/1908.03831 https://inspirehep.net/search?p=find+EPRINT+arXiv:1908.03831 https://doi.org/10.1007/JHEP08(2018)003 https://doi.org/10.1007/JHEP08(2018)003 https://arxiv.org/abs/1803.11089 https://inspirehep.net/search?p=find+EPRINT+arXiv:1803.11089 https://doi.org/10.1016/j.physletb.2018.04.046 https://arxiv.org/abs/1802.02445 https://inspirehep.net/search?p=find+EPRINT+arXiv:1802.02445 https://doi.org/10.1140/epjc/s10052-018-5557-y https://arxiv.org/abs/1706.01473 https://inspirehep.net/search?p=find+EPRINT+arXiv:1706.01473 https://doi.org/10.1007/JHEP06(2019)028 https://doi.org/10.1007/JHEP06(2019)028 https://arxiv.org/abs/1902.08474 https://inspirehep.net/search?p=find+EPRINT+arXiv:1902.08474 https://doi.org/10.1016/j.physletb.2018.11.037 https://arxiv.org/abs/1805.00736 https://inspirehep.net/search?p=find+EPRINT+arXiv:1805.00736 https://doi.org/10.1103/PhysRevD.93.011502 https://arxiv.org/abs/1509.06392 https://inspirehep.net/search?p=find+EPRINT+arXiv:1509.06392 https://doi.org/10.1103/PhysRevD.93.054004 https://arxiv.org/abs/1511.05590 https://inspirehep.net/search?p=find+EPRINT+arXiv:1511.05590 https://doi.org/10.1007/JHEP07(2018)172 https://arxiv.org/abs/1805.07243 https://inspirehep.net/search?p=find+EPRINT+arXiv:1805.07243 https://doi.org/10.1016/j.nuclphysb.2005.12.022 https://arxiv.org/abs/hep-ph/0508068 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0508068 https://doi.org/10.1103/PhysRevD.81.093007 https://arxiv.org/abs/0911.4135 https://inspirehep.net/search?p=find+EPRINT+arXiv:0911.4135 https://doi.org/10.1007/JHEP11(2011)064 https://arxiv.org/abs/1109.2109 https://inspirehep.net/search?p=find+EPRINT+arXiv:1109.2109 https://doi.org/10.1007/JHEP12(2018)132 https://arxiv.org/abs/1805.05916 https://inspirehep.net/search?p=find+EPRINT+arXiv:1805.05916 https://doi.org/10.1103/PhysRevLett.106.132001 https://arxiv.org/abs/1011.4225 https://inspirehep.net/search?p=find+EPRINT+arXiv:1011.4225 https://doi.org/10.1103/PhysRevD.84.051503 https://doi.org/10.1103/PhysRevD.84.051503 https://arxiv.org/abs/1105.1991 https://inspirehep.net/search?p=find+EPRINT+arXiv:1105.1991 https://doi.org/10.1103/PhysRevD.85.045003 https://arxiv.org/abs/1109.6293 https://inspirehep.net/search?p=find+EPRINT+arXiv:1109.6293 J H E P 1 1 ( 2 0 1 9 ) 1 2 1 [34] C. Pisano, D. Boer, S.J. Brodsky, M.G.A. Buffing and P.J. Mulders, Linear polarization of gluons and photons in unpolarized collider experiments, JHEP 10 (2013) 024 [arXiv:1307.3417] [INSPIRE]. [35] A. Dumitru, T. Lappi and V. Skokov, Distribution of Linearly Polarized Gluons and Elliptic Azimuthal Anisotropy in Deep Inelastic Scattering Dijet Production at High Energy, Phys. Rev. Lett. 115 (2015) 252301 [arXiv:1508.04438] [INSPIRE]. [36] A. Mukherjee and S. Rajesh, Linearly polarized gluons in charmonium and bottomonium production in color octet model, Phys. Rev. D 95 (2017) 034039 [arXiv:1611.05974] [INSPIRE]. [37] D. Boer, P.J. Mulders, J. Zhou and Y.-j. Zhou, Suppression of maximal linear gluon polarization in angular asymmetries, JHEP 10 (2017) 196 [arXiv:1702.08195] [INSPIRE]. [38] A.V. Efremov, N.Y. Ivanov and O.V. Teryaev, How to measure the linear polarization of gluons in unpolarized proton using the heavy-quark pair leptoproduction, Phys. Lett. B 777 (2018) 435 [arXiv:1711.05221] [INSPIRE]. [39] A.V. Efremov, N.Y. Ivanov and O.V. Teryaev, The ratio R = dσL/dσT in heavy-quark pair leptoproduction as a probe of linearly polarized gluons in unpolarized proton, Phys. Lett. B 780 (2018) 303 [arXiv:1801.03398] [INSPIRE]. [40] R. Kishore and A. Mukherjee, Accessing linearly polarized gluon distribution in J/ψ production at the electron-ion collider, Phys. Rev. D 99 (2019) 054012 [arXiv:1811.07495] [INSPIRE]. [41] M.G. Echevarria, Proper TMD factorization for quarkonia production: pp→ ηc,b as a study case, JHEP 10 (2019) 144 [arXiv:1907.06494] [INSPIRE]. [42] C. Marquet, C. Roiesnel and P. Taels, Linearly polarized small-x gluons in forward heavy-quark pair production, Phys. Rev. D 97 (2018) 014004 [arXiv:1710.05698] [INSPIRE]. [43] L.D. McLerran and R. Venugopalan, Gluon distribution functions for very large nuclei at small transverse momentum, Phys. Rev. D 49 (1994) 3352 [hep-ph/9311205] [INSPIRE]. [44] artemide web-page, https://teorica.fis.ucm.es/artemide/. [45] artemide repository, https://github.com/vladimirovalexey/artemide-public. [46] I. Scimemi, A. Tarasov and A. Vladimirov, Collinear matching for Sivers function at next-to-leading order, JHEP 05 (2019) 125 [arXiv:1901.04519] [INSPIRE]. [47] T. Gehrmann, T. Lubbert and L.L. Yang, Transverse parton distribution functions at next-to-next-to-leading order: the quark-to-quark case, Phys. Rev. Lett. 109 (2012) 242003 [arXiv:1209.0682] [INSPIRE]. [48] M.G. Echevarŕıa, A. Idilbi and I. Scimemi, Soft and Collinear Factorization and Transverse Momentum Dependent Parton Distribution Functions, Phys. Lett. B 726 (2013) 795 [arXiv:1211.1947] [INSPIRE]. [49] S. Moch, J.A.M. Vermaseren and A. Vogt, The Quark form-factor at higher orders, JHEP 08 (2005) 049 [hep-ph/0507039] [INSPIRE]. [50] T. Gehrmann, E.W.N. Glover, T. Huber, N. Ikizlerli and C. Studerus, Calculation of the quark and gluon form factors to three loops in QCD, JHEP 06 (2010) 094 [arXiv:1004.3653] [INSPIRE]. – 22 – https://doi.org/10.1007/JHEP10(2013)024 https://arxiv.org/abs/1307.3417 https://inspirehep.net/search?p=find+EPRINT+arXiv:1307.3417 https://doi.org/10.1103/PhysRevLett.115.252301 https://doi.org/10.1103/PhysRevLett.115.252301 https://arxiv.org/abs/1508.04438 https://inspirehep.net/search?p=find+EPRINT+arXiv:1508.04438 https://doi.org/10.1103/PhysRevD.95.034039 https://arxiv.org/abs/1611.05974 https://inspirehep.net/search?p=find+EPRINT+arXiv:1611.05974 https://doi.org/10.1007/JHEP10(2017)196 https://arxiv.org/abs/1702.08195 https://inspirehep.net/search?p=find+EPRINT+arXiv:1702.08195 https://doi.org/10.1016/j.physletb.2017.12.061 https://doi.org/10.1016/j.physletb.2017.12.061 https://arxiv.org/abs/1711.05221 https://inspirehep.net/search?p=find+EPRINT+arXiv:1711.05221 https://doi.org/10.1016/j.physletb.2018.03.008 https://doi.org/10.1016/j.physletb.2018.03.008 https://arxiv.org/abs/1801.03398 https://inspirehep.net/search?p=find+EPRINT+arXiv:1801.03398 https://doi.org/10.1103/PhysRevD.99.054012 https://arxiv.org/abs/1811.07495 https://inspirehep.net/search?p=find+EPRINT+arXiv:1811.07495 https://doi.org/10.1007/JHEP10(2019)144 https://arxiv.org/abs/1907.06494 https://inspirehep.net/search?p=find+EPRINT+arXiv:1907.06494 https://doi.org/10.1103/PhysRevD.97.014004 https://arxiv.org/abs/1710.05698 https://inspirehep.net/search?p=find+EPRINT+arXiv:1710.05698 https://doi.org/10.1103/PhysRevD.49.3352 https://arxiv.org/abs/hep-ph/9311205 https://inspirehep.net/search?p=find+EPRINT+hep-ph/9311205 https://teorica.fis.ucm.es/artemide/ https://github.com/vladimirovalexey/artemide-public https://doi.org/10.1007/JHEP05(2019)125 https://arxiv.org/abs/1901.04519 https://inspirehep.net/search?p=find+EPRINT+arXiv:1901.04519 https://doi.org/10.1103/PhysRevLett.109.242003 https://arxiv.org/abs/1209.0682 https://inspirehep.net/search?p=find+EPRINT+arXiv:1209.0682 https://doi.org/10.1016/j.physletb.2013.09.003 https://arxiv.org/abs/1211.1947 https://inspirehep.net/search?p=find+EPRINT+arXiv:1211.1947 https://doi.org/10.1088/1126-6708/2005/08/049 https://doi.org/10.1088/1126-6708/2005/08/049 https://arxiv.org/abs/hep-ph/0507039 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0507039 https://doi.org/10.1007/JHEP06(2010)094 https://arxiv.org/abs/1004.3653 https://inspirehep.net/search?p=find+EPRINT+arXiv:1004.3653 J H E P 1 1 ( 2 0 1 9 ) 1 2 1 [51] A.A. Vladimirov, Correspondence between Soft and Rapidity Anomalous Dimensions, Phys. Rev. Lett. 118 (2017) 062001 [arXiv:1610.05791] [INSPIRE]. [52] Y. Li and H.X. Zhu, Bootstrapping Rapidity Anomalous Dimensions for Transverse-Momentum Resummation, Phys. Rev. Lett. 118 (2017) 022004 [arXiv:1604.01404] [INSPIRE]. [53] J.C. Collins and D.E. Soper, Back-To-Back Jets: Fourier Transform from B to K-Transverse, Nucl. Phys. B 197 (1982) 446 [INSPIRE]. [54] I. Scimemi and A. Vladimirov, Power corrections and renormalons in Transverse Momentum Distributions, JHEP 03 (2017) 002 [arXiv:1609.06047] [INSPIRE]. [55] M.-X. Luo, T.-Z. Yang, H.X. Zhu and Y.J. Zhu, Transverse Parton Distribution and Fragmentation Functions at NNLO: the Gluon Case, arXiv:1909.13820 [INSPIRE]. [56] Y. Li, D. Neill and H.X. Zhu, An Exponential Regulator for Rapidity Divergences, arXiv:1604.00392 [INSPIRE]. [57] M.A. Shifman, A.I. Vainshtein, M.B. Voloshin and V.I. Zakharov, Low-Energy Theorems for Higgs Boson Couplings to Photons, Sov. J. Nucl. Phys. 30 (1979) 711 [INSPIRE]. [58] M. Krämer, E. Laenen and M. Spira, Soft gluon radiation in Higgs boson production at the LHC, Nucl. Phys. B 511 (1998) 523 [hep-ph/9611272] [INSPIRE]. [59] K.G. Chetyrkin, B.A. Kniehl and M. Steinhauser, Hadronic Higgs decay to order α4 S , Phys. Rev. Lett. 79 (1997) 353 [hep-ph/9705240] [INSPIRE]. [60] H.M. Georgi, S.L. Glashow, M.E. Machacek and D.V. Nanopoulos, Higgs Bosons from Two Gluon Annihilation in Proton Proton Collisions, Phys. Rev. Lett. 40 (1978) 692 [INSPIRE]. [61] V. Ahrens, T. Becher, M. Neubert and L.L. Yang, Renormalization-Group Improved Prediction for Higgs Production at Hadron Colliders, Eur. Phys. J. C 62 (2009) 333 [arXiv:0809.4283] [INSPIRE]. [62] V. Ravindran, J. Smith and W.L. van Neerven, NNLO corrections to the total cross-section for Higgs boson production in hadron hadron collisions, Nucl. Phys. B 665 (2003) 325 [hep-ph/0302135] [INSPIRE]. [63] C. Anastasiou, K. Melnikov and F. Petriello, Fully differential Higgs boson production and the di-photon signal through next-to-next-to-leading order, Nucl. Phys. B 724 (2005) 197 [hep-ph/0501130] [INSPIRE]. [64] V. Ahrens, T. Becher, M. Neubert and L.L. Yang, Origin of the Large Perturbative Corrections to Higgs Production at Hadron Colliders, Phys. Rev. D 79 (2009) 033013 [arXiv:0808.3008] [INSPIRE]. [65] NNPDF collaboration, Parton distributions from high-precision collider data, Eur. Phys. J. C 77 (2017) 663 [arXiv:1706.00428] [INSPIRE]. [66] A. Vladimirov, Pion-induced Drell-Yan processes within TMD factorization, JHEP 10 (2019) 090 [arXiv:1907.10356] [INSPIRE]. [67] T. Sjöstrand, S. Mrenna and P.Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05 (2006) 026 [hep-ph/0603175] [INSPIRE]. [68] T. Sjöstrand, S. Mrenna and P.Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput. Phys. Commun. 178 (2008) 852 [arXiv:0710.3820] [INSPIRE]. – 23 – https://doi.org/10.1103/PhysRevLett.118.062001 https://doi.org/10.1103/PhysRevLett.118.062001 https://arxiv.org/abs/1610.05791 https://inspirehep.net/search?p=find+EPRINT+arXiv:1610.05791 https://doi.org/10.1103/PhysRevLett.118.022004 https://arxiv.org/abs/1604.01404 https://inspirehep.net/search?p=find+EPRINT+arXiv:1604.01404 https://doi.org/10.1016/0550-3213(82)90453-9 https://inspirehep.net/search?p=find+J+%22Nucl.Phys.,B197,446%22 https://doi.org/10.1007/JHEP03(2017)002 https://arxiv.org/abs/1609.06047 https://inspirehep.net/search?p=find+EPRINT+arXiv:1609.06047 https://arxiv.org/abs/1909.13820 https://inspirehep.net/search?p=find+EPRINT+arXiv:1909.13820 https://arxiv.org/abs/1604.00392 https://inspirehep.net/search?p=find+EPRINT+arXiv:1604.00392 https://inspirehep.net/search?p=find+J+%22Sov.J.Nucl.Phys.,30,711%22 https://doi.org/10.1016/S0550-3213(97)00679-2 https://arxiv.org/abs/hep-ph/9611272 https://inspirehep.net/search?p=find+EPRINT+hep-ph/9611272 https://doi.org/10.1103/PhysRevLett.79.353 https://doi.org/10.1103/PhysRevLett.79.353 https://arxiv.org/abs/hep-ph/9705240 https://inspirehep.net/search?p=find+EPRINT+hep-ph/9705240 https://doi.org/10.1103/PhysRevLett.40.692 https://inspirehep.net/search?p=find+J+%22Phys.Rev.Lett.,40,692%22 https://doi.org/10.1140/epjc/s10052-009-1030-2 https://arxiv.org/abs/0809.4283 https://inspirehep.net/search?p=find+EPRINT+arXiv:0809.4283 https://doi.org/10.1016/S0550-3213(03)00457-7 https://arxiv.org/abs/hep-ph/0302135 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0302135 https://doi.org/10.1016/j.nuclphysb.2005.06.036 https://arxiv.org/abs/hep-ph/0501130 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0501130 https://doi.org/10.1103/PhysRevD.79.033013 https://arxiv.org/abs/0808.3008 https://inspirehep.net/search?p=find+EPRINT+arXiv:0808.3008 https://doi.org/10.1140/epjc/s10052-017-5199-5 https://doi.org/10.1140/epjc/s10052-017-5199-5 https://arxiv.org/abs/1706.00428 https://inspirehep.net/search?p=find+EPRINT+arXiv:1706.00428 https://doi.org/10.1007/JHEP10(2019)090 https://doi.org/10.1007/JHEP10(2019)090 https://arxiv.org/abs/1907.10356 https://inspirehep.net/search?p=find+EPRINT+arXiv:1907.10356 https://doi.org/10.1088/1126-6708/2006/05/026 https://doi.org/10.1088/1126-6708/2006/05/026 https://arxiv.org/abs/hep-ph/0603175 https://inspirehep.net/search?p=find+EPRINT+hep-ph/0603175 https://doi.org/10.1016/j.cpc.2008.01.036 https://doi.org/10.1016/j.cpc.2008.01.036 https://arxiv.org/abs/0710.3820 https://inspirehep.net/search?p=find+EPRINT+arXiv:0710.3820 J H E P 1 1 ( 2 0 1 9 ) 1 2 1 [69] CMS collaboration, Measurement of differential cross sections for Higgs boson production in the diphoton decay channel in pp collisions at √ s = 8 TeV, Eur. Phys. J. C 76 (2016) 13 [arXiv:1508.07819] [INSPIRE]. [70] V.A. Smirnov, Evaluating Feynman integrals, Springer Tracts Mod. Phys. 211 (2004) 1. [71] G. Altarelli and G. Parisi, Asymptotic Freedom in Parton Language, Nucl. Phys. B 126 (1977) 298 [INSPIRE]. – 24 – https://doi.org/10.1140/epjc/s10052-015-3853-3 https://arxiv.org/abs/1508.07819 https://inspirehep.net/search?p=find+EPRINT+arXiv:1508.07819 https://doi.org/10.1007/b95498 https://doi.org/10.1016/0550-3213(77)90384-4 https://doi.org/10.1016/0550-3213(77)90384-4 https://inspirehep.net/search?p=find+J+%22Nucl.Phys.,B126,298%22 Introduction Gluon TMD distributions Definition OPE at small-b Renormalization of TMDPDF Matching coefficient for lpTMDPDF at NNLO Evaluation of the matching coefficient Logarithmic part of the coefficient function Finite part of the coefficient function lpTMDPDF at NNLO and its contribution Higgs production Conclusions Relevant set of master integrals for linearly polarized gluon TMD Logarithm terms of matching coefficient for lpTMDPDF