Universality of transverse-momentum resummation and hard factors at the NNLO

We consider QCD radiative corrections to the production of colourless high-mass systems in hadron collisions. The logarithmically-enhanced contributions at small transverse momentum are treated to all perturbative orders by a universal resummation formula that depends on a single process-dependent hard factor. We show that the hard factor is directly related to the all-order virtual amplitude of the corresponding partonic process. The direct relation is universal (process independent), and it is expressed by an all-order factorization formula that we explicitly evaluate up to the next-to-next-to-leading order (NNLO) in QCD perturbation theory. Once the NNLO scattering amplitude is available, the corresponding hard factor is directly determined: it controls NNLO contributions in resummed calculations at full next-to-next-to-leading logarithmic accuracy, and it can be used in applications of the q_T subtraction formalism to perform fully-exclusive perturbative calculations up to NNLO. The universality structure of the hard factor and its explicit NNLO form are also extended to the related formalism of threshold resummation.


Introduction
The transverse-momentum (q T ) distribution of systems with high invariant mass M produced in hadron collisions is important for physics studies within and beyond the Standard Model (SM). This paper is devoted to a theoretical study of QCD radiative corrections to transverse-momentum distributions.
We consider the inclusive hard-scattering reaction where the collision of the two hadrons h 1 and h 2 with momenta p 1 and p 2 produces the triggered final state F , and X denotes the accompanying final-state radiation. The observed final state F is a generic system of one or more colourless particles, such as lepton pairs (produced by Drell-Yan (DY) mechanism), photon pairs, vector bosons, Higgs boson(s), and so forth. The momenta of these final state particles are denoted by q 1 ,q 2 ...q n . The system F has total invariant mass M 2 = (q 1 + q 2 + ...q n ) 2 , transverse momentum q T and rapidity y. We use √ s to denote the centre-of-mass energy of the colliding hadrons, which are treated in the massless approximation (s = (p 1 + p 2 ) 2 = 2p 1 · p 2 ).
The transverse-momentum cross section for the process in Eq. (1) is computable by using perturbative QCD. However, in the small-q T region (roughly, in the region where q T ≪ M) the convergence of the fixed-order perturbative expansion in powers of the QCD coupling α S is spoiled by the presence of large logarithmic terms of the type ln n (M 2 /q 2 T ). The predictivity of perturbative QCD can be recovered through the summation of these logarithmically-enhanced contributions to all order in α S [1,2,3,4,5,6,7,8,9,10]. As already stated, we shall limit ourselves to considering the production of systems F of non-strongly interacting particles. The all-order analysis of the q T distribution of systems F that involve coloured QCD partons has just started to be investigated, by considering [11] the specific case in which F is formed by a tt pair.
In the case of a generic system F of colourless particles, the large logarithmic contributions to the q T cross section can be systematically resummed to all perturbative orders, and the structure of the resummed calculation can be organized in a process-independent form [4,6,9,10]. The allorder resummation formalism was first developed for the DY process [6] (and the kinematicallyrelated process of two-particle correlations in e + e − annihilation [4]). The process-independent extension of the formalism has required two additional main steps: the understanding of the allorder process-independent structure of the Sudakov form factor (through the factorization of a single process-dependent hard factor) [9], and the complete generalization to processes that are initiated by the gluon fusion mechanism [10].
The all-order process-independent form of the resummed calculation has a factorized structure, whose resummation factors are (see Sect. 2) the (quark and gluon) Sudakov form factor, processindependent collinear factors and a process-dependent hard or, more precisely (see Sect. 4), hardvirtual factor. The resummation of the logarithmic contributions is controlled by these factors or, equivalently, by a corresponding set of perturbative functions whose perturbative resummation coefficients are computable order-by-order in α S . The perturbative coefficients of the Sudakov form factor are known, since some time [5,7,12,8,13], up to the second order in α S , and the third-order coefficient A (3) (which is necessary to explicitly perform resummation up to the nextto-next-to-leading logarithmic (NNLL) accuracy) is also known [14]. The next-to-next-to-leading order (NNLO) QCD calculation of the q T cross section (in the small-q T region) has been explicitly carried out in analytic form for two benchmark processes, namely, SM Higgs boson production [15] and the DY process [16]. The results of Refs. [15,16] provide us with the complete knowledge of the process-independent collinear resummation coefficients up to the second order in α S , and with the explicit expression of the hard coefficients for these two specific processes. The purpose of the present paper (see below) is to explicitly point out and derive the underlying universal (processindependent) structure of the process-dependent hard factor of the QCD all-order resummation formalism.
In Refs. [17,18,14,19,20,21,22], the resummation of small-q T logarithms has been reformulated in terms of factorization formulae that involve Soft Collinear Effective Theory operators and (process-dependent) hard matching coefficients. The formulation of Ref. [14] has been applied [23] to the DY process by explicitly computing the (process-independent) collinear quark-quark coefficients and the DY hard coefficient at the NNLO. The results of this calculation [23] agree with those obtained in Ref. [16]. Transverse-momentum cross sections can also be studied by using other approaches (which go beyond the customary QCD resummation formalism of the present paper) that use transverse-momentum dependent (TMD) factorization (see Refs. [24,25,26,27] and references therein) and, consequently, k ⊥ -unintegrated parton densities and partonic cross sections that are both TMD quantities.
In this paper we study the process-dependent hard factor of the transverse-momentum resummation formula. We show that, for any process of the class in Eq. (1), the all-order hard factor has a universal structure that involves a minimal amount of process-dependent information. The process-dependent information is entirely given by the scattering amplitude of the Born-level partonic subprocess and its virtual radiative corrections. Knowing the scattering amplitude, the hard-virtual resummation factor is determined by a universal (process-independent) factorization formula. The universality structure of the factorization formula has a soft (and collinear) origin, and it is closely (though indirectly) related to the universal structure of the infrared divergences [28] of the scattering amplitude. This process-independent structure of the hard-virtual term, which generalizes the next-to-leading order (NLO) results of Ref. [13], is valid to all perturbative orders, and we explicitly determine the process-independent form of the hard-virtual term up to the NNLO. Using this general NNLO result, the hard-virtual resummation factor for each process of the class in Eq. (1) is straightforwardly computable up to its NNLO, provided the corresponding scattering amplitude is known.
In the final part of the paper, we consider the related formalism of threshold resummation [29,30] for the total cross section. The process-independent formalism of threshold resummation also involves a corresponding process-dependent hard factor. We shall show that this factor has a universality structure that is analogous to the case of transverse-momentum resummation. In particular, we directly relate the process-dependent hard factors for transverse-momentum and threshold resummation in a form that is fully universal and completely independent of each specific process (e.g., independent of the corresponding scattering amplitude).
The knowledge of the NNLO hard-virtual term completes the q T resummation formalism in explicit form up to full NNLL+NNLO accuracy. This permits direct applications to NNLL+NNLO resummed calculations for any processes of the class in Eq. (1) (provided the corresponding NNLO amplitude is known), as already done for the cases of SM Higgs boson [31,32,33] and DY [34,35] production.
The NNLO information of the q T resummation formalism is also relevant in the context of fixed order calculations. Indeed, it permits to carry out fully-exclusive NNLO calculations by applying the q T subtraction formalism of Ref. [36] (the subtraction counterterms of the formalism follow [36] from the fixed-order expansion of the q T resummation formula, as in Sect. 2.4 of Ref. [31]). The q T subtraction formalism has been applied to the NNLO computation of Higgs boson [36,37] and vector boson production [38], associated production of the Higgs boson with a W boson [39], diphoton production [40] and Zγ production [41]. The computations of Refs. [36,37,38,39] were based on the specific calculation of the NNLO hard-virtual coefficients of the corresponding processes [15,16]. The computations of Refs. [40,41] used the NNLO hard-virtual coefficients that are determined by applying the universal form of the hard-virtual term that is derived and illustrated in the present paper.
The paper is organized as follows. In Sect. 2 we recall the transverse-momentum resummation formalism in impact parameter space, and we introduce our notation. In Sect. 3 we present the explicit expressions of the process-independent resummation coefficients up to NNLO. Section 4 is devoted to the process-dependent hard coefficients. We discuss and illustrate the universal allorder form of the hard-virtual coefficients by relating them to the process-dependent scattering amplitudes, through the introduction of suitably subtracted hard-virtual matrix elements. The process-independent structure of the hard-virtual coefficients is explicitly computed up to the NNLO. In Sect. 5 we extend our discussion and results on the universal structure of the hardvirtual coefficients to the case of threshold resummation. In Sect. 6 we summarize our results. In Appendix A we report the explicit expressions of the NLO and NNLO hard-virtual coefficients for DY, Higgs boson and diphoton production.
We consider the inclusive-production process in Eq. (1), and we introduce the corresponding fully differential cross section which depends on the total momentum of the system F (i.e. on the variables q T , M, y). The cross section also depends on the set of additional variables that controls the kinematics of the particles in the system F . In Eq. (2) these additional variables are generically denoted as Ω = {Ω A , Ω B , . . . } (correspondingly, we define dΩ ≡ dΩ A dΩ B . . . ). To be general, we do not explicitly specify these variables, and we only require that the kinematical variables {Ω A , Ω B , . . . } are independent of q T , M and y and that the set of variables {q T , M, y, Ω} completely determines the kinematical configuration (i.e., the momenta q i ) of the particles in the system F . For instance, if the system F is formed by two particles, there are only two variables in the set Ω, and they can be the rapidity y i and the azimuthal angle φ(q T i ) of one of the two particles (the particle with momentum q i ). Note that the cross section in Eq. (2) and the corresponding resummation formula can be straightforwardly integrated with respect to some of the final-state variables {Ω A , Ω B , . . . }, thus leading to results for observables that are more inclusive than the differential cross section in Eq. (2).
We also recall that we are considering the production of a system F of colourless particles or, more precisely, a system of non-strongly interacting particles (i.e., F cannot include QCD partons and their fragmentation products). Therefore, at the Born (lowest-order) level, the cross section in Eq. (2) is controlled by the partonic subprocesses of quark-antiquark (qq) annihilation, and gluon fusion, Owing to colour conservation, no other partonic subprocesses can occur at the Born level. More importantly (see below), the distinction between qq annihilation and gluon fusion leads to relevant (and physical) differences [42,10] in the context of small-q T resummation.
To study the q T dependence of the differential cross section in Eq. (2) within QCD perturbation theory, we introduce the following decomposition: Both terms in the right-hand side are obtained through convolutions of partonic cross sections and the scale-dependent parton distributions f a/h (x, µ 2 ) (a = q f ,q f , g is the parton label) of the colliding hadrons. We use parton densities as defined in the MS factorization scheme, and α S (q 2 ) is the QCD running coupling in the MS renormalization scheme. The partonic cross sections that enter the singular component (the first term in the right-hand side of Eq. (5)) contain all the contributions that are enhanced (or 'singular') at small q T . These contributions are proportional to δ (2) (q T ) or to large logarithms of the type 1 q 2 T ln m (M 2 /q 2 T ). On the contrary, the partonic cross sections of the second term in the right-hand side of Eq. (5) are regular (i.e. free of logarithmic terms) order-by-order in perturbation theory as q T → 0. More precisely, the integration of dσ (reg) F /d 2 q T over the range 0 ≤ q T ≤ Q 0 leads to a finite result that, at each fixed order in α S , vanishes in the limit Q 0 → 0.
The regular component dσ (reg) F of the q T cross section depends on the specific process in Eq. (1) that we are considering. In the following we focus on the singular component, dσ (sing) F , which has a universal all-order structure. The corresponding resummation formula is written as [6,9,10] where b 0 = 2e −γ E (γ E = 0.5772 . . . is the Euler number) is a numerical coefficient, and the kinematical variables x 1 and x 2 are The right-hand side of Eq. (6) involves the Fourier transformation with respect to the impact parameter b and two convolutions over the longitudinal-momentum fractions z 1 and z 2 . The parton densities f a i /h i (x, µ 2 ) of the colliding hadrons are evaluated at the scale µ = b 0 /b, which depends on the impact parameter. The function S c (M, b) is the Sudakov form factor. This factor, which only depends on the type (c = q or c = g) of colliding partons, is universal (process independent) [9], and it resums the logarithmically-enhanced contributions of the form ln M 2 b 2 (the region q T ≪ M corresponds to Mb ≫ 1 in impact parameter space). The all-order expression of S c (M, b) is [6] [14] are explicitly known.
The factor that is symbolically denoted by dσ (0) cc, F in Eq. (6) is the Born-level cross section dσ (0) (i.e., the cross section at its corresponding lowest order in α S ) of the partonic subprocesses cc → F in Eqs. (3) and (4) (in the case of the cc = qq annihilation channel, the quark and antiquark can actually have different flavours). Making the symbolic notation explicit, we have where x 1 p µ 1 (x 2 p µ 2 ) is the momentum of the parton c (c). In Eq. (6), we have included the contribution of both the qq annihilation channel (c = q,q) and the gluon fusion channel (c = g); one of these two contributing channels may be absent (i.e. dσ (0) cc, F = 0 in that channel), depending on the specific final-state system F .
The Born level factor dσ (0) cc, F is obviously process dependent, although its process dependence is elementary (it is simply due to the Born level scattering amplitude of the partonic process cc → F ). The remaining process dependence of Eq. (6) is embodied in the 'hard-collinear' factor H F C 1 C 2 . This factor includes a process-independent part and a process-dependent part. The structure of the process-dependent part is the main subject of the present paper.
In the case of processes that are initiated at the Born level by the qq annihilation channel (c = q), the symbolic factor H F C 1 C 2 in Eq. (6) has the following explicit form [9] and the functions H F q and C q a = Cqā have the perturbative expansion The function H F q is process dependent, whereas the functions C q a are universal (they only depend on the parton indices). We add an important remark on the factorized structure in the right-hand side of Eq. (11): the scale of α S is M 2 in the case of H F q , whereas the scale is b 2 0 /b 2 in the case of C q a . The appearance of these two different scales is essential [9] to disentangle the process dependence of H F q from the process-independent Sudakov form factor (S q ) and collinear functions (C q a ).
In the case of processes that are initiated at the Born level by the gluon fusion channel (c = g), the physics of the small-q cross section has a richer structure, which is the consequence of collinear correlations [10] that are produced by the evolution of the colliding hadrons into gluon partonic states. In this case, the resummation formula (6) and, specifically, its factor H F C 1 C 2 are more involved than those for the qq channel, since collinear radiation from the colliding gluons leads to spin and azimuthal correlations [42,10]. The symbolic factor H F C 1 C 2 in Eq. (6) has the following explicit form [10]: where the function H F g has the perturbative expansion and the following lowest-order normalization: Analogously to Eq. (11), in Eq. (14) the function H F g;µ 1 ν 1 ,µ 2 ν 2 is process dependent (and it is controlled by α S at the scale M 2 ) and the partonic functions C µ ν ga are process independent (and they are controlled by α S at the scale b 2 0 /b 2 ). At variance with Eq. (11) (where the factorization structure in the right-hand side is independent of the degrees of freedom of the colliding quark and antiquarks), in Eq. (14) the process-dependent function H F g depends on the Lorentz indices (and, thus, on the spins) {µ i ν i } of the colliding gluons with momenta x i p i (i = 1, 2) and this dependence is coupled to (and correlated with) a corresponding dependence of the partonic functions C µ i ν i ga i . The Lorentz tensor coefficients C µ i ν i ga i in Eq. (14) depend on b 2 (through the scale of α S ) and, moreover, they also depend on the direction (i.e., the azimuthal angle) of the impact parameter vector b in the transverse plane. The structure of the partonic tensor is [10] C µν g a (z; where and b µ = (0, b, 0) is the two-dimensional impact parameter vector in the four-dimensional notation . The gluonic coefficient function C g a (z; α S ) (a = q,q, g) in the right-hand side of Eq. (17) has the same perturbative structure as in Eq. (13), and it reads In contrast, the perturbative expansion of the coefficient functions G ga , which are specific to gluon-initiated processes, starts at O(α S ), and we write We recall [10] an important physical consequence of the different small-q T resummation structure between the qq annihilation and gluon fusion channels: the absence of azimuthal correlations with respect to q T in the qq annihilation channel, and the presence of correlations with a definite predictable azimuthal dependence in the gluon fusion channel. Indeed, in the case of qq annihilation, all the factors in the integrand of the Fourier transformation on the right-hand side of the resummation formula (6) are functions of b 2 , with no dependence on the azimuthal angle φ(b) of b. Therefore, the integration over φ(b) in Eq. (6) can be straightforwardly carried out, and it leads [4,6] to a one-dimensional Bessel transformation that involves the 0th-order Bessel function J 0 (bq T ). This implies that the right-hand side of Eq. (6) and, hence, the singular part of the q T differential cross section depend only on q 2 T , with no additional dependence on the azimuthal angle φ(q T ) of q T . Unlike the case of qq annihilation, the gluon fusion factor H F C 1 C 2 in Eqs. (6) and (14) does depend on the azimuthal angle φ(b) of the impact parameter b. Therefore, the integration over φ(b) in the Fourier transformation of Eq. (6) is more complicated. It leads to one-dimensional Bessel transformations that involve J 0 (bq T ) and higher-order Bessel functions, such as the 2-nd order and 4th-order functions J 2 (bq T ) and J 4 (bq T ). More importantly, it leads to a definite structure of azimuthal correlations with respect to the azimuthal angle φ(q T ) of the transverse momentum q T . The small-q T cross section in Eq. (6) can be expressed [10] in terms of a contribution that does not depend on φ(q T ) plus a contribution that is given by a linear combination of the four angular functions cos (2φ(q T )) , sin (2φ(q T )) , cos (4φ(q T )) and sin (4φ(q T )). No other functional dependence on φ(q T ) is allowed by the resummation formula (6) in the gluon fusion channel.
We recall that, due to its specific factorization structure, the resummation formula in Eq. (6) is invariant under the following renormalization-group transformation [9] where h c (α S ) = 1 + O(α S ) is an arbitrary perturbative function (with h q (α S ) = hq(α S )). More precisely, in the case of gluon-initiated processes, Eq. (24) becomes In the right-hand side of Eq. (23), β(α S ) denotes the QCD β-function: where N f is the number of quark flavours, N c is the number of colours, and the colour factors are As a consequence of the renormalizationgroup symmetry in Eqs. (22)-(25), the resummation factors H F , S c , C qa and C µν ga are not separately defined (and, thus, computable) in an unambiguous way. Equivalently, each of these separate factors can be precisely defined only by specifying a resummation scheme [9].
To present the main results of this paper in the following sections, we find it convenient to specify a resummation scheme. Therefore, in the rest of this paper we work in the scheme, dubbed hard scheme, that is defined as follows. The flavour off-diagonal coefficients C gg (z) in Eqs. (13) and (20) is instead due to both 'regular' functions and 'singular' distributions in the limit z → 1. The 'singular' distributions are δ(1 − z) and the customary plus-distributions of the form [(ln k (1 − z))/(1 − z)] + (k = 0, 1, 2 . . . ). The hard scheme is the scheme in which, order-by-order in perturbation theory, the coefficients C (n) ab (z) with n ≥ 1 do not contain any δ(1 − z) term. We remark (see also Sect. 4) that this definition directly implies that all the process-dependent virtual corrections to the Born level subprocesses in Eqs. (3) and (4) are embodied in the resummation coefficient H F c .
We note that the specification of the hard scheme (or any other scheme) has sole practical purposes of presentation (theoretical results can be equivalently presented, as actually done in Refs. [15] and [16], by explicitly parametrizing the resummation-scheme dependence of the resummation factors). Having presented explicit results in the hard scheme, they can be translated in other schemes by properly choosing the functions h c (α S ) (c = q, g) and applying the transformation in Eqs. (22)- (25). Moreover, and more importantly, the q T cross section, its all-order resummation formula (6) and any consistent perturbative truncation (either order-by-order in α S or in classes of logarithmic terms) of the latter [9,31] are completely independent of the resummation scheme.
The process-independent partonic coefficients C ab (z; α S ) in Eqs. (11) and (14) are explicitly known up to the NNLO (see references in Sect. 3). The universality structure of the processdependent coefficients H F c at NNLO and higher orders (see Sect. 4) is one of the main result of the present work.

Process-independent coefficients
Before discussing the general structure of the resummation coefficients H F c , in this section we present the expressions of the process-independent resummation coefficients in the hard scheme, which is defined in Sec. 2.
The partonic functions C ab and G gb in Eqs. (13), (20) and (21) depend on the parton indices. Owing to charge conjugation invariance and flavour symmetry of QCD, the dependence on the parton indeces is fully specified by the five independent quark functions {C qq , C qq ′ , C qq , C qq ′ , C qg } [16] (q and q ′ denote quarks with different flavour) and the four independent gluon functions The first-order coefficients C (1) ab (z) are explicitly known [12,43,8,13]. Their expressions in the hard scheme can be obtained from their corresponding expression in an arbitrary scheme by simply setting the coefficient of the δ(1 − z) term to zero. We get The first-order coefficients G (1) ga are resummation-scheme independent, and they read [10] where C a is the Casimir colour coefficient of the parton a with C q = C F and C g = C A .
The second-order process-independent collinear coefficients C (2) ab (z) of Eqs. (13) and (20) have been computed in Refs. [36,38,15,16]. The quark-quark coefficient C (2) qq (z) has been independently computed in Ref. [23]. The expressions of these coefficients in the hard scheme can be straightforwardly obtained from the results of Refs. [15,16] and are explicitly reported below.
The second-order gluon collinear coefficients G (2) ga (z) (a = q, g) of Eq. (21) are not yet known. We can comment on the role of G (2) ga in practical terms. In the specific and important case of Higgs boson production by gluon fusion, the coefficient G ga does not contribute to the cross section at the NNLO (and NNLL accuracy). The Higgs boson cross section is discussed in detail in Ref. [10]: by direct inspection of Eq. (45) of Ref. [10], we can see that G (2) ga starts to contribute at the N 3 LO. In most of the other processes (e.g., F = γγ, Zγ, , W + W − ), the system F can be produced by both qq annihilation and gluon fusion. In these case, due to the absence of direct coupling of the gluons to the colourless particles in the system F , the production channel gg → F is suppressed by some powers of α S with respect to the channel qq → F . Therefore, also in these cases the coefficient G ga does not contribute to the NNLO cross section. This formal conclusion (based on counting the powers of α S ) has a caveat, since the gg → F channel can receive a quantitative enhancement from the possibly large luminosity of the gluon parton densities. However, the knowledge of the first-order coefficients C (1) ga and G (1) ga should be sufficient to compute the contribution from the channel gg → F to a quantitative level that is comparable to that of the contribution from the channel qq → F (whose collinear coefficients are fully known up to the second order). In summary, we conclude that the effect of G (2) ga rarely contributes in actual (practical) computations of the q T cross section at the NNLO or NNLL accuracy.

Hard-virtual coefficients
In this section we focus on the process-dependent coefficient H F . In the hard scheme that we are using, this coefficient contains all the information on the process-dependent virtual corrections, and, therefore, we can show that H F can be related in a process-independent (universal) way to the multiloop virtual amplitude M cc→F of the partonic process cc → F . In the following we first specify the notation that we use to denote the all-loop virtual amplitude M cc→F . Then we introduce an auxiliary (hard-virtual) amplitude M cc→F that is directly obtained from M cc→F by using a process-independent relation. Finally, we use the hard-virtual amplitude M cc→F to present the explicit expression of the hard-virtual coefficient H F up to the NNLO.
We consider the partonic elastic-production process where the two colliding partons with momentap 1 andp 2 are either cc = gg or cc = qq (we do not explicitly denote the flavour of the quark q, although in the case with cc = qq, the quark and the antiquark can have different flavours), and F ({q i }) is the triggered final-state system in Eq. (1). The loop scattering amplitude of the process in Eq. (41) contains ultraviolet (UV) and infrared (IR) singularities, which are regularized in d = 4 − 2ǫ space-time dimensions by using the customary scheme of conventional dimensional regularization (CDR) † . Before performing renormalization, the multiloop QCD amplitude has a perturbative dependence on powers of α u S µ 2ǫ 0 , where α u S is the bare coupling and µ 0 is the dimensional-regularization scale. In the following we work with the renormalized on-shell scattering amplitude that is obtained from the corresponding unrenormalized amplitude by just expressing the bare coupling α u S in terms of the running coupling α S (µ 2 R ) according to the MS scheme relation where µ R is the renormalization scale, β 0 and β 1 are the first two coefficients of the QCD β-function in Eq. (28) and the factor S ǫ is The renormalized all-loop amplitude of the process in Eq. (41) is denoted by M cc→F (p 1 ,p 2 ; {q i }), and it has the perturbative (loop) expansion where the value k of the overall power of α S depends on the specific process (for instance, k = 0 in the case of the vector boson production process qq → V , and k = 1 in the case of the Higgs boson production process gg → H through a heavy-quark loop). Note also that the lowest-order perturbative term M (0) cc→F is not necessarily a tree-level amplitude (for instance, it involves a quark loop in the cases gg → H and gg → γγ). The perturbative terms M (l) cc→F (l = 1, 2, . . . ) are UV finite, but they still depend on ǫ (although this dependence is not explicitly denoted in Eq. (44)). In particular, the amplitude M (l) cc→F at the l-th perturbative order is IR divergent as ǫ → 0, and it behaves as where the dots stand for ǫ-poles of lower order. The IR divergent contributions to the scattering amplitude have a universal structure [28], which is explicitly known at the one-loop [47,28], two-loop [28,48] and three-loop [49,50] level for the class of processes in Eq. (41).
The explicit calculations and the results of Ref. [13] show that the NLO hard-virtual coefficient H F (1) is explicitly related in a process-independent form to the leading-order (LO) amplitude M (0) cc→F and to the IR finite part of the NLO amplitude M (1) cc→F . The relation between H F c and M cc→F can be extended to the NNLO and to higher-order levels. This extension can be formulated and expressed in simple and general terms by introducing an auxiliary (hard-virtual) amplitude M cc→F that is directly obtained from M cc→F in a universal (process-independent) way. In practice, M cc→F is obtained from M cc→F by removing its IR divergences and a definite amount of IR finite terms. The (IR divergent and finite) terms that are removed from M cc→F originate from real emission contributions to the cross section and, therefore, these terms and M cc→F specifically depend on the transverse-momentum cross section of Eq. (2), which we consider throughout this paper.
The hard-virtual amplitude M cc→F has a perturbative expansion that is analogous to that in Eq. (44). We write At the LO, M cc→F and M cc→F coincide, and we have At higher-perturbative orders, M In Eqs. (48) and (49) c act as IR subtraction operators (factors), and their functional dependence is explicitly denoted in Eqs. (48) and (49). These terms are process independent (they do not depend on F and on its specific production mechanism in Eq. (41)): they only depend on the invariant mass M of the system F (through the dimensionless ratio M 2 /µ 2 R ), on the type c (c = q, g) of colliding partons, and on ǫ. In particular, I cc→F are IR finite as ǫ → 0. We also note that the structure of Eqs. (48) and (49) and the explicit dependence ofĨ The explicit expression of the first-order (one-loop) subtraction operatorĨ (1) a is and The coefficient δ q T affects only the IR finite part of the subtraction operator. The known results on the NLO hard-collinear coefficients H F (1) c [13] are recovered by fixing The second-order (two-loop) subtraction operatorĨ with where H (2) coll a is the contribution that is proportional to γ a (1) and H (2) soft a is the remaining contribution (which is proportional to C a ) in Eq. (57). The QCD coefficients K in Eq. (55) and d (1) in Eq. (57) (they control the IR divergences ofĨ (2) a ) are [28] and the coefficients γ a (1) (a = q, g) are given in Eqs. (35) and (36). The coefficient δ q T (1) in Eq. (57) affects only the IR finite part of the two-loop subtraction operator. We find (see Sect. 4.1) Having introduced the subtracted amplitude M cc→F , we can relate it to the process-dependent resummation coefficients H F c of Eqs. (6), (11) and (14). In the case of processes initiated by qq annihilation (see Eqs. (11) and (12)), the all-order coefficient H F q can be written as where k is the value of the overall power of α S in the expansion of M cc→F (see Eqs. (44) and (46)). In the case of processes initiated by gluon fusion (see Eqs. (14) and (15)), the analogue of Eq. (61) is and the all-order coefficient H F g is [10] where d µν = d µν (p 1 , p 2 ) is the polarization tensor in Eq. (18) and it projects onto the Lorentz indices in the transverse plane.
In Eqs. (61) and (62), the notation | M qq→F | 2 and |M  do not depend on the spin. We recall that, according to the notation in Eq. (2), the kinematics of the final-state momenta {q i } is fully specified by the total momentum q = i q i and the set of variables Ω. Therefore, the dependence of the amplitudes on {q i } completely determines the Ω dependence of H F c in Eqs. (61) and (62). To be precise, we also note that H F c is computed in d = 4 space-time dimensions and, therefore, the right-hand side of Eqs. (61) and (62) has to be evaluated in the limit ǫ → 0 (this limit is well-defined and straightforward, since M (0) cc→F and the order-by-order expansion of the hard-virtual amplitude M cc→F are IR finite). An additional remark regards the dependence on the renormalization scale. According to the resummation formula (6), the all-order factor dσ (0) cc, F H F c (and, consequently, the left-hand side of Eqs. (61) and (62)) is renormalization-group invariant, and it is perturbatively computable as series expansion in powers of α S (M 2 ), with no dependence on µ R . This property is fully consistent with the form of Eqs. (61) and (62), since the all-order hard-virtual amplitude in the right-hand side of these equations is a renormalization-group invariant quantity. Obviously, each side of these equations can be expanded in powers of α S (µ 2 R ), thus leading to corresponding perturbative coefficients that explicitly depend on M 2 /µ 2 R .

The structure of the hard-virtual term
The all-loop amplitude M cc→F receives contributions from loop momenta in different kinematical regions. Roughly speaking, the loop momentum k can be in the UV region (k ≫ M), in the IR region (k ≪ M) or in the hard (intermediate) region (k ∼ M). The UV contributions are treated and 'removed' by renormalization. As already anticipated (and discussed below), the subtraction factorĨ c has a IR (soft and collinear) origin, and it 'removes' the IR contributions from M cc→F . The subtracted (hard-virtual) amplitude M cc→F can thus be interpreted as originating from the hard component of the virtual radiative corrections to the LO amplitude M with The factorization formula (64) gives the all-order relation between the hard-virtual amplitude M cc→F (which determines H F c through Eqs. (61) and (62)) and the scattering amplitude M cc→F (the perturbative expansion of Eqs. (64) up to the NNLO exactly gives Eqs. (47)-(49)). The allorder subtraction factorĨ c (ǫ, M 2 ) in Eq. (64) is independent of µ R and, thus, it is renormalizationgroup invariant. The order-by-order dependence on µ R simply arises from the expansion (see Eq. (65)) in terms of powers of α S (µ 2 R ) and perturbative coefficientsĨ (n) c (ǫ, M 2 /µ 2 R ). Note that I c (ǫ, M 2 ) depends on ǫ, and we are referring to renormalization-group invariance in d = 4 − 2ǫ dimensions (i.e., to all orders in ǫ), where the right-hand side of Eq. (26) has to be modified by replacing β(α S ) with the d-dimensional β-function β(ǫ, α S ) = − ǫ + β(α S ). The µ R independence ofĨ c (ǫ, M 2 ) can be explicitly checked up to the NNLO by using the expressions ofĨ  We can illustrate the origin and the derivation of the results in Eqs. (61), (62) and (64) by starting from the direct computation of the q T cross section in Eq. (2). The calculation of the cross section or, more precisely, of the corresponding partonic cross sections involves three types of contributions: (i) the elastic-production process in Eq. (41); (ii) inelastic (real-emission) processes, where the system F is accompanied by additional final-state partons; (iii) the collinear counterterm that is necessary to define the MS parton densities in terms of the bare (naïve) parton densities.
(i) The elastic process directly contributes to the partonic cross section and thus to H F c with a term that is proportional to (the square of) the all-loop amplitude M cc→F .
(ii) Since we are interested in the small-q T singular cross section of Eq. (6), the calculation of the inelastic processes can be simplified. In the small-q T limit, the additional final-state partons in the inelastic processes must be either soft or collinear to one of the colliding partons (non-soft and non-collinear partons give cross section contributions that are relatively suppressed by some powers of q T /M ∼ 1/(bM)). The radiation of soft [51,52,53] and collinear [54,51,52,55] partons from two colliding partons is described by QCD factorization formulae, where the singular soft/collinear term (which includes its virtual radiative corrections) is universal (process independent) and it acts onto the all-loop amplitude M cc→F as in the factorized expression on the right-hand side of Eq. (64). To be precise, soft/collinear factorization works at the amplitude (and squared amplitude) level, and it can be spoiled by kinematical effects at the cross section level, i.e., after the (cross section dependent) phase-space integration of the squared amplitudes. However, in our case, factorization breaking effects of kinematical origin cannot arise, since we are effectively working in impact parameter space (in the small-q T limit, the kinematics of the q T cross section is exactly factorized [2] by the Fourier transformation to b space). More precisely, we are considering the cross section contributions at fixed values of the impact parameter b, namely, the integrand terms on the right-hand side of Eq. (6): the factorized structure of these terms directly follows form soft/collinear factorization formulae. The following step in our discussion consists in the observation that the radiation of collinear or, more precisely, non-soft collinear partons requires a non-vanishing longitudinal-momentum recoil and, therefore, it cannot contribute to the factor H F c in Eq. (6) (non-soft collinear radiation definitely contribute to the other b-dependent factors on the right-hand side of Eq. (6)). In summary, considering the inelastic processes, the factor H F c receives contributions only from soft radiation: these are the factorized soft contributions that are b-independent and that do not vanish in the near-elastic limit.
(iii) The calculation of the elastic and inelastic processes gives the bare partonic cross section. The introduction of the (IR divergent) collinear counterterm of the parton densities amounts to multiply the entire bare partonic cross section with a process-independent factor that has a convolution structure with respect to the longitudinal-momentum fractions of the colliding partons. The evolution kernel of the convolution is the Altarelli-Parisi splitting function P ab (z; α S ), and only the soft part (analogously to the contribution of the inelastic processes) and the virtual part (i.e., the part that is proportional to δ(1 − z)) of the splitting function can contribute to the factor H F c in Eq. (6). This 'virtual-collinear' contribution to H F c is process independent (it only depends on the type, i.e. quarks or gluons, of colliding partons), and it has the same factorized structure as in Eq. (64). Note that the explicit form of the collinear counterterm is fully specified [44], since we are considering parton densities defined in the MS factorization scheme. Moreover, we can add that the entire virtual part of the collinear counterterm is included in H F c , since we are working in the hard scheme (where the perturbative corrections to the coefficients C ab (z; α S ) of Eq. (6) contain no contributions that are proportional to δ(1 − z)).
In summary, from our general discussion we can conclude that, in the hard scheme, the resummation factor H F c has the structure given by Eqs. (61), (62) and (64). The all-loop amplitude M cc→F of Eq. (64) originates from the elastic process in Eq. (41). The remaining universal factor I c (ǫ, M 2 ) on the right-hand side of Eq. (64) includes two types of contributions: a soft contribution (from inelastic processes and the collinear counterterm) and a collinear contribution from the virtual part of the MS collinear counterterm. In the following, we combine these conclusions with two additional properties of H F c (its renormalization-group invariance and its IR finiteness), and we shall show that the explicit form ofĨ c (ǫ, M 2 ) is (almost) completely determined up to the NNLO.
The IR divergences of the all-loop amplitude M cc→F have a known universal structure that can be presented in the following form [28]: where the all-loop factor I c has a perturbative expansion that is analogous to that in Eq. (65). The component M fin. cc→F of the amplitude is IR finite as ǫ → 0, while the process-independent factor I c includes IR divergent ǫ-poles. The perturbative (loop) expansion of Eq. (66) iteratively determines the IR divergent component of M cc→F (see Eqs. (12) and (18) in Ref. [28]). Note that the separation between M fin. cc→F and IR divergent terms depends on the amount of IR finite contributions that are actually included in I c (in particular, as discussed below, the hard-virtual amplitude M cc→F can be regarded as a specific definition of M fin. cc→F ).
The relation (66) can be rewritten as and it can be directly compared with the form of the hard-virtual amplitude M cc→F . The relations (64) and (67) are in one-to-one correspondence through the simple replacement M cc→F ↔ M fin.
cc→F andĨ c ↔ I c . Therefore, by requiring that M cc→F is IR finite, we conclude that the ǫ-pole contributions of the operatorĨ c (ǫ, M 2 ) are exactly the same as those of I c (equivalently,Ĩ c and I c can differ only through terms that produce IR finite contributions to M cc→F ).
We continue our all-order discussion at the fixed-order level to make our conclusions more explicit and clear. We consider the NLO and NNLO terms. The extension to higher orders is straightforward.
At the NLO (one-loop) level,Ĩ    48)). Then, we note that the entire expression on the right-hand side of Eq. (52) exactly coincides with the virtual part of the collinear counterterm in the MS factorization scheme [44] (we recall that in the MS factorization scheme the collinear counterterm has only ǫ-pole contributions, and that the coefficient γ a in Eqs. (52) and (53) is equal to the coefficient of the virtual part of the first-order Altarelli-Parisi splitting function). Therefore, we can conclude that the only unknown contribution toĨ (1) a is an ǫ-independent term that has a soft origin. This term is included in the right-hand side of Eq. (51), and it can be written as δ a = C a δ q T , since the intensity of soft radiation from the parton a is simply proportional to the Casimir coefficient C a of that parton. In summary, we have given a proof of the results in Eqs. (50)-(53), although we cannot give the explicit value of the process-independent coefficient δ q T on the basis of our general discussion. The explicit determination of δ q T requires a detailed calculation, and the value in Eq. (54) is taken from available results in the literature [13].
At the NNLO (two-loop) level we can repeat the same reasoning and steps as at the NLO (one-loop) level. Using the explicit form of the ǫ poles of the operator I (2) a [28] in Eq. (67) and the requirement of renormalization-group invariance, we eventually obtain the second-order subtraction operatorĨ (2) coll a is completely determined apart from ǫ-independent contributions of soft and collinear (from the virtual part of the collinear counterterm) origin. However, the ǫ-independent contribution of collinear origin is vanishing (as in the case ofĨ coincides with the entire collinear-counterterm contribution (in the MS factorization scheme) due to the virtual part of the second-order Altarelli-Parisi splitting function [44,45]. Therefore, the remaining ǫ-independent contribution to H (2) a has a soft origin, and it can be written in the form δ a(1) = C a δ q T (1) in the right-hand side of Eq. (57). Owing to its origin from soft factorization (see the Appendix of Ref. [51] and Sect. 5 of Ref. [53]), δ a(1) is simply proportional to the Casimir coefficient C a of the radiating (colliding) parton a and the QCD coefficient δ q T (1) is fully process independent, namely, it is the same coefficient for processes that are initiated by either qq annihilation or gluon fusion.
In summary, we have proven the two-loop results in Eqs. (55)-(59), although δ q T (1) cannot be determined from our general discussion. The explicit determination of δ q T (1) requires a detailed calculation. Such a calculation can be explicitly performed in a general process-independent form by extending the analysis of Ref. [13] (which is based on NNLO soft/collinear factorization formulae [51,52,53,54,55]) to the necessary level of accuracy (in practice, by including contributions at higher order in ǫ that were omitted in the actual computation of Ref. [13]). Alternatively, we can exploit our proof of the universality of δ q T (1) and, therefore, we can determine the value of δ q T from the NNLO calculation of a single specific process. We have followed the latter procedure (see below) to obtain the explicit result of δ q T (1) that is reported in Eq. (60). The NNLO computation of the DY cross section at small values of q T was performed in Ref. [38], and the complete result is presented in Ref. [16] in explicit analytic form. From Ref. [16] we can thus extract the explicit value of the NNLO coefficient H DY (2) q for the DY process (see Appendix A). The same value of H DY (2) q is obtained by the fully independent calculation of Ref. [23]. The scattering amplitude M qq→DY for the DY process was computed long ago up to the two-loop level [56,57]. Therefore, using the result of Refs. [56,57]  (1) as an unknown parameter, we thus obtain an expression of H DY (2) q (which linearly depends on δ q T (1) ) that can be directly compared with the explicit value extracted from the calculation of Refs. [16,23]. This comparison gives the value of δ q T (1) that is reported in Eq. (60).
The same procedure can be applied to extract the value of δ q T (1) from Higgs boson production by gluon fusion. Indeed, also for this process we know both the explicit value of the coefficient H H (2) g from a direct NNLO computation [15] (see Appendix A) and the corresponding two-loop amplitude M gg→H [48] (both results use the large-M top approximation). Using these results, we confirm the value of δ q T (1) that we have extracted from the DY process.
Note that the agreement between these two independent determinations (extractions) of δ q T is a highly non-trivial check of the results of Sect. 4, especially because we are considering two processes that are controlled by the qq annihilation channel and the gluon fusion channel (δ q T (1) is instead independent of the specific channel). Note also that this agreement can alternatively (i.e., assuming the knowledge of δ q T (1) ) be regarded as a non-trivial (though partial) cross-check of the results of the NNLO calculations of DY [16,23] and Higgs boson [15] production.
We add a final comment related to our general discussion on the structure of the hard-virtual term. The all-loop scattering amplitude M cc→F includes an overall phase factor e +iφ Coul. (ǫ,M 2 ) (the phase φ Coul. (ǫ, M 2 ) is IR divergent), which is the QCD analogue of the QED Coulomb phase. This phase factor is physically harmless, since it cancels in the evaluation of the squared amplitude and, consequently, in the computation of cross sections. Our explicit expression ofĨ (1) c includes an imaginary contribution (the term that is proportional to iπ/ǫ in Eq. (51)), and corresponding contributions are present in the expression (55) ofĨ (2) c (through its dependence onĨ (1) c ). These contributions of 'imaginary' origin exactly correspond to the perturbative expansion of the Coulomb phase factor, and they lead to a hard-virtual amplitude M cc→F that does not include the overall and harmless (but IR divergent) factor e +iφ Coul. (ǫ,M 2 ) . This imaginary contributions toĨ c cannot arise from a direct computation of H F c at the cross section level and, actually, we have introduced them to the sole practical (aesthetical) purpose of cancelling the IR divergent Coulomb phase of M cc→F . In other words, by removing these contributions fromĨ

Universality and threshold resummation
The structure of transverse-momentum resummation and, especially, of the hard-virtual term can be compared with the analogous structure of threshold resummation [29,30], which arises in the context of the QCD computation of the total cross section. To highlight the main aspects of the comparison, we consider the total cross section for the process of Eq. (1) in the simple case (the restriction to this simple case has the sole purpose of simplifying the notation) in which the final-state system F consists of a single ('on-shell') particle of mass M (for example, F can be a vector boson or a Higgs boson). The total cross section σ F (p 1 , p 2 ; M 2 ) for the production of the system F is computable in QCD perturbation theory according to the following factorization formula: whereσ F a 1 a 2 is the total partonic cross section for the inclusive partonic process a 1 a 2 → F + X and, for simplicity, the parton densities f a i /h i (z i , M 2 ) (i = 1, 2) are evaluated at the scale M 2 (the inclusion of an arbitrary factorization scale µ F in the parton densities and in the partonic cross sections can be implemented in a straightforward way by using the Altarelli-Parisi evolution equations of f a/h (z, µ 2 F )). The partonic cross sectionσ F a 1 a 2 (ŝ; M 2 ; α S (M 2 )) depends on the mass M of the system F , on the centre-of-mass energyŝ of the colliding partons, and it is a renormalizationgroup invariant quantity that can be perturbatively computed as series expansion in powers of α S (M 2 ) (equivalently, we can expandσ F a 1 a 2 in powers of α S (µ 2 R ), with corresponding perturbative coefficients that explicitly depend on M 2 /µ 2 R ).
The kinematical ratio z = M 2 /ŝ parametrizes the distance from the partonic threshold. In the kinematical region close to the partonic threshold (i.e., where z → 1), the partonic cross sectionσ F a 1 a 2 receives large QCD radiative corrections of the type 1 1−z ln m (1 − z) + (the subscript '+' denotes the customary 'plus-distribution'). The all-order resummation of these logarithmic contributions can be systematically performed by working in Mellin (N-moment) space [29,30].
The Mellin transformσ N (M 2 ) of the partonic cross sectionσ(ŝ; M 2 ) is defined aŝ In Mellin space, the threshold region z → 1 corresponds to the limit N → ∞, and the plusdistributions become powers of ln N ( 1 1−z ln m (1 − z) + → ln m+1 N + 'subleading logs ′ ). These logarithmic contributions are evaluated to all perturbative orders by using threshold resummation [29,30]. Neglecting terms that are relatively suppressed by powers of 1/N in the limit N → ∞, we writeσ cc, N has a universal (process-independent) all-order structure that is given by the following threshold-resummation formula [29,30,58,59]: where σ (0) cc→F is the lowest-order cross section for the partonic process cc → F . The radiative factor ∆ c, N resums all the perturbative contributions α n S ln m N (including some constant terms, i.e. terms with m = 0). This factor only depends on the type (c = q or c = g) of colliding partons (∆ c, N does not depend on the final-state system F ), and it has the form where A th c (α S ) and D c (α S ) are perturbative series in α S , The perturbative coefficients A  The factor C th cc→F in Eq. (71) embodies remaining N-independent contributions (i.e., terms that are constant in the limit N → ∞) to the partonic cross section. This factor is definitely process dependent, and it has the general perturbative expansion The NLO and NNLO coefficients C th (1) and C th (2) are explicitly known in the case of DY [57,59] and Higgs boson [67,64,68,58] production. Considering these two specific processes, relations between the N-independent factor C th (α S ) and the corresponding virtual amplitudes (namely, the quark and gluon form factors) were discussed and examined in Refs. [69,65,66,70,71]. Considering a generic process of the class in Eq. (1) and using soft-gluon factorization formulae [51,53,28], the authors of Ref. [72] have recently shown how the NLO and NNLO coefficients C th (1) cc→F and C th (2) cc→F can be directly and explicitly related in a process-independent form to the oneloop and two-loop scattering amplitude M cc→F of the underlying partonic process in Eq. (41).
The results in Eqs. (75)-(79) are obtained by using the same reasoning and discussion as in Sect. 4.1 and, in particular, by exploiting the properties of soft/collinear factorization. We do not repeat the entire discussion of Sect. 4.1, and we limit ourselves to remarking the few points in which the discussion slightly differs. Considering the computation of the total partonic cross section, we can directly refer to the classification in contributions from (i) the elastic-production process, (ii) inelastic processes and (iii) the collinear counterterm. (i) The elastic process directly contributes to M th cc→F in Eq. (76) (C th cc→F in Eq. (75)) with the all-loop scattering amplitude M cc→F (with the squared amplitude). (ii) In the kinematical region close to the partonic threshold, the inelastic processes contribute to the partonic cross section of Eq. (71) only through final-state radiation of soft partons (in the case of transverse-momentum resummation, collinear radiation also contributes, and it is responsible for the presence of the collinear coefficients C qa and C µν ga in Eq. (6) and for the differences between the qq and gg channels). Soft factorization at the (squared) amplitude level is not spoiled by kinematical effects, since the kinematics of the total cross section is exactly factorized [29,30] by the Mellin transformation to N space. This leads to the same conclusion as in Sect. 4.1 about the contribution of the inelastic processes to the hard-virtual term. This contribution is factorized and it has a soft origin in both cases of transverse-momentum and threshold resummation. In particular, at the cross section level (i.e., after the corresponding phase space integration), this soft term produces contributions to the coefficients δ th and δ th (1) in the expressions ofĨ th c and C th cc→F , analogously to the corresponding contributions to the coefficients δ q T and δ q T (1) in the expressions ofĨ c and H F c . (iii) The radiative factor ∆ c, N in Eq. (71) is entirely due to soft radiation [29,30]. Therefore, the complete virtual part of the collinear counterterm in the MS factorization scheme directly contributes to M th cc→F and C th cc→F (analogously to the contribution to M cc→F and H F c ). It follows that the collinear counterterm contributions toĨ th c andĨ c are completely analogous.
Owing to their process independence, the threshold resummation coefficients δ th and δ th (1) (analogously to δ q T and δ q T (1) ) can be explicitly evaluated from either the NNLO calculation of a single process or an NNLO calculation in a process-independent form. We have explicitly verified that the general results in Eqs. (75)- (77) and the explicit values of the coefficients in Eqs. (78) and (79) are consistent with the NNLO results of the process-independent calculation of Ref. [72].
The close correspondence between the hard-virtual terms H F c and C th cc→F of the resummation formulae in Eqs. (6) and (71) can also be expressed in direct form. Using Eqs. (61)-(63), Eq. (75) and the expressions (64) and (76) of the corresponding hard-virtual amplitudes, we obtain where, in the case of gluon fusion processes, the numerator in the left-hand side of Eq. (80) is defined as H F g ≡ g µ 1 ν 1 g µ 2 ν 2 H F µ 1 ν 1 ,µ 2 ν 2 g . The equality in Eq. (81) is obtained by using the explicit expression ofĨ c (ǫ, M 2 ) (see Eqs. (50)- (59)) and the corresponding expression ofĨ th c (ǫ, M 2 ) (see Eq. (77) and accompanying comments) up to the NNLO.
Considering the ratio of hard-virtual terms for a specific process as in Eq. (80), the effect of the all-loop amplitude M cc→F (and the associated process dependence) entirely cancels. This ratio is completely determined by the contribution of the inelastic processes (namely, the factorized radiation of final-state partons and the corresponding virtual corrections) to the corresponding cross sections. The ensuing IR (soft and collinear) singularities completely cancel, and the final expression in Eq. (81) is entirely determined by the IR finite contributions due to the (real) emission of soft QCD radiation. The exponent in Eq. (81) is directly proportional to the Casimir factor C c (i.e., the colour charge) of the colliding partons c andc : this proportionality is a straightforward consequence of the exponentiating correlation structure of the factorization formulae for softparton radiation from QCD squared amplitudes [51,53]. The value of the perturbative coefficients δ q T − δ th and δ q T (1) − δ th (1) in the exponent has a kinematical origin.

Summary
In this paper we have considered QCD radiative corrections to the production of a generic colourless high-mass system F in hadronic collisions (see Eq. (1)). Large logarithmic terms arise in the QCD perturbative expansion when the high-mass system F is produced at small transverse momentum. These logarithmic terms can be resummed to all perturbative orders by using a universal (process-independent) resummation formula (see Sect. 2) and, then, they are controlled by a set of resummation factors and ensuing perturbative resummation coefficients. After having recalled the process independence of the Sudakov form factor and the explicit expressions (up to NNLO) of the process-independent collinear coefficients (see Sect. 3), in Sect. 4 we have focused on the hard-virtual factor H F c . We have shown that, although this factor is process dependent, it can be directly related (see Eqs. (61)-(63)) in a universal (process-independent) way to the IR finite part of the all-order virtual amplitude M cc→F of the corresponding partonic subprocess cc → F . Therefore, the all-order scattering amplitude M cc→F is the sole process-dependent information that is eventually required by the all-order resummation formula. The relation between H F c and M cc→F follows from a universal all-order factorization formula (see Eqs. (64) and (65)) that originates from the factorization properties of soft (and collinear) parton radiation. We have explicitly determined this relation up to the NNLO. More precisely, we have shown that this relation is fully determined by the structure of IR singularities of the all-order amplitude M cc→F and by renormalization-group invariance up to a single coefficient (of soft origin) at each perturbative order. We have explicitly determined these coefficients at NLO and NNLO. Therefore, knowing the NNLO scattering amplitude M cc→F , its corresponding hard-virtual resummation factor H F c is straightforwardly determined up to NNLO.
The results presented in this paper, with the knowledge of the other process-independent resummation coefficients (which are recalled in Sects. 2 and 3), complete (modulo the second-order coefficients G (2) ga , as discussed at the end of Sect. 3) the q T resummation formalism in explicit form up to full NNLL and NNLO accuracy for all the processes in the class of Eq. (1). This permits applications to NNLL+NNLO resummed calculations for any processes whose NNLO scattering amplitudes are available. Moreover, since the hard-virtual (and collinear) resummation coefficients are exactly the coefficients that are required to implement the q T subtraction formalism [36], the results that we have presented are directly and straightforwardly applicable to perform fullyexclusive NNLO computations for each of these processes.
We have also considered the related process-independent formalism of threshold resummation (see Sect. 5). We have shown that the process-dependent hard-virtual factor C th cc→F of threshold resummation has a universal (process-independent) structure that is analogous to that of the hard-virtual factor H F c of transverse-momentum resummation. The process-independent relation between C th cc→F and the scattering amplitude M cc→F has been explicitly pointed out up to NNLO. In particular, we have shown that, for each specific process, the ratio H F c /C th cc→F is completely independent of the process (i.e., independent of M cc→F ), and it is fully determined by the associated soft-parton radiation.

A Hard-virtual coefficients in DY, Higgs boson and diphoton production
In this Appendix we report the explicit expressions of the hard-virtual coefficients in the hard scheme for the cases of DY, Higgs boson and diphoton production. The NNLO coefficients H DY (2) q and H H(2) g for DY and Higgs boson (using the large-m Q approximation) production were obtained in Refs. [16] and [15] by performing a direct QCD computation of the corresponding q T cross sections. The same results can be recovered (as discussed in Sect. 4.1) by using the processindependent structure of the hard-virtual coefficients. The explicit expression (which is presented in Eq. (90)) of the NNLO coefficient H γγ(2) q for diphoton production is directly obtained by using the results of Sect. 4.
Starting from the DY process, we consider the production (through the qq annihilation channel) of a virtual photon or a vector boson (V = γ * , W ± , Z) and the subsequent leptonic decay. The corresponding Born-level cross section dσ in the hard scheme are [16] H DY (1) We then consider the production of the SM Higgs boson H through the gluon fusion channel, where H couples to a heavy-quark loop. If the Higgs boson decays non-hadronically (e.g., H → ZZ → 4 leptons, or H → γγ), the dependence on the kinematics of its decay products only affects the Born-level cross section dσ (0) gg,H in Eq. (6), while the corresponding hard-virtual factor only depends on α S (M 2 ) (as in the case of DY production) and on the mass of the heavy quark in the loop. The same conclusion applies to the hadronic decay of H, if we neglect the QCD interferences between the initial and final states. In both cases (i.e., non-hadronic decay or hadronic decay without interferences), the spin (Lorentz index) correlation structure of the hard-collinear factor in Eq. (14) can be simplified. Indeed, as shown in Ref. [10], the right-hand side of Eq. (14) turns out to be proportional to the following (Lorentz scalar) hard-virtual factor: H H g (α S (M 2 )) = g µ 1 ν 1 g µ 2 ν 2 H H µ 1 ν 1 ,µ 2 ν 2 g (α S (M 2 )) , whose perturbative coefficients H H (n) g follow from the perturbative expansion in Eq. (15). Assuming that the Higgs boson couples to a single heavy quark of mass m Q , the first-order coefficient H The expression in Eq. (85) is obtained by using the process-independent formulae in Eqs. (48), (62) and (63) In the large-m Q limit, the explicit expression of the NNLO hard-virtual coefficient H H(2) g in the hard scheme is [15] H H ( where L Q = ln(M 2 /m 2 Q ). The scattering amplitude M gg→H has been computed [76] up to its NNLO by including corrections to the large-m Q approximation (the evaluation of the corrections uses the expansion parameter 1/m 2 Q ). Using the process-independent formulae in Eqs. (49), (62) and (63), these corrections can be straightforwardly included in the expression of the NNLO hard-virtual coefficient H We finally consider the case of diphoton production. In this case, the production cross section receives contributions from both the qq annihilation and gluon fusion channels. The final-state system F = γγ can be produced by the elastic partonic subprocesses qq → γγ and gg → γγ. In the perturbative evaluation of the cross section, the subprocess qq → γγ first contributes at the LO (through the corresponding tree-level scattering amplitude), while the subprocess gg → γγ starts to contributes at the NNLO (through the corresponding scattering amplitude that involves the one-loop QCD interaction of light and heavy quarks). Therefore, to the purpose of evaluating the complete q T cross section up to its NNLO, the gluon fusion channel only contributes at its corresponding lowest order. Having computed dσ (0) gg,γγ , the NNLO contribution to Eq. (6) from the gluon fusion channel is obtained by simply considering the lowest-order hard-virtual coefficient H γγ(0) g (in practice, we can simply implement the replacement H γγ µ 1 ν 1 ,µ 2 ν 2 g → H γγ(0)µ 1 ν 1 ,µ 2 ν 2 g → d µ 1 ν 1 (p 1 , p 2 ) d µ 2 ν 2 (p 1 , p 2 )/4 in Eqs. (14) and (15)). The hard-virtual coefficient H γγ q of the subprocess qq → γγ has instead to be explicitly evaluated up to its NNLO (i.e., we need the perturbative coefficients H γγ(1) q and H γγ(2) q ).
Using the notation of Eq. (41), to compute H γγ q we have to consider the partonic process q(p 1 )q(p 2 ) → γ(q 1 )γ(q 2 ), whose Mandelstam kinematical variables arê s = (p 1 +p 2 ) 2 = M 2 ,û = (p 2 − q 1 ) 2 ,t = (p 1 − q 2 ) 2 , with the constraintŝ +t +û = M 2 +t +û = 0. Unlike the case of DY and Higgs boson production, the hard-virtual term H γγ q of Eq. (11) depends on the kinematical variables of the final-state diphoton system. We explicitly specify these two variables (they are generically denoted as Ω = {Ω A , Ω B } in Eqs. (2), (6) and (11)) by using the azimuthal angle of q 1 and the ratio v = −û/ŝ = −û/M 2 . The hard-virtual term H γγ q is invariant with respect to azimuthal rotations (see Eq. (61)) and it is a dimensionless function of its kinematical variables. This implies that H γγ q only depends on v, and we simply use the notation H γγ q (v; α S (M 2 )). The first-order coefficient H γγ(1) q (v) of Eq. (12) is known [77]. Its explicit expression in the hard scheme is The second-order coefficient H γγ(2) q (v) was computed and used in Ref. [40]. The calculation of H γγ(2) q is performed by using the universality structure of the hard-virtual term (see Eqs. (47)- (49) and (61)) and the explicit result [78] of the two-loop amplitude M qq→γγ of the process qq → γγ ‡ . The result that we obtain for H γγ (2) where the functions F 0×2 inite,qqγγ;s and F 1×1 inite,qqγγ;s are defined in Eqs. (4.6) and (5.3) of Ref. [78], respectively, and the function A LO (v) is (91) ‡ We note that there are some typos in Ref. [78]: in Eq. (3.13) the factor Γ(1 − ǫ)/Γ(1 − 2ǫ) has to be replaced with the factor Γ(1 − 2ǫ)/Γ(1 − ǫ); the overall sign on the right-hand side of Eqs. (C.1), (C.2) and (C.3) has to be reversed (in particular, without this flip of sign, the coefficients of the IR poles of the NLO and NNLO scattering amplitude have wrong signs).