Transverse-momentum resummation and the spectrum of the Higgs boson at the LHC

We consider the transverse-momentum (q_T) distribution of generic high-mass systems (lepton pairs, vector bosons, Higgs particles, ....) produced in hadron collisions. At small q_T, we concentrate on the all-order resummation of the logarithmically-enhanced contributions in QCD perturbation theory. We elaborate on the $b$-space resummation formalism and introduce some novel features: the large logarithmic contributions are systematically exponentiated in a process-independent form and, after integration over q_T, they are constrained by perturbative unitarity to give a vanishing contribution to the total cross section. At intermediate and large q_T, resummation is consistently combined with fixed-order perturbative results, to obtain predictions with uniform theoretical accuracy over the entire range of transverse momenta. The formalism is applied to Standard Model Higgs boson production at LHC energies. We combine the most advanced perturbative information available at present for this process: resummation up to next-to-next-to-leading logarithmic accuracy and fixed-order perturbation theory up to next-to-leading order. The results show a high stability with respect to perturbative QCD uncertainties.


Introduction
to reduce the associated perturbative uncertainty at the level of about ±10% [17].
When studying the q T distribution of the Higgs boson in QCD perturbation theory, it is convenient to start by considering separately the large-q T and small-q T regions.
The large-q T region is identified by the condition q T ∼ M H . In this region, the perturbative series is controlled by a small expansion parameter, α S (M 2 H ), and calculations based on the truncation of the series at a fixed order in α S are theoretically justified. SM Higgs boson production at large q T via gluon fusion has to be accompanied by the radiation of at least one recoiling parton, so the LO term for this observable is of O(α 3 S ). The LO calculation was reported in Ref. [18]; it shows that the large-M t approximation works well as long as M H ∼ < 2M t and q T ∼ < M t . Similar results on the validity of the large-M t approximation were obtained in the case of the associated production of a Higgs boson plus 2 jets (2 recoiling partons at large transverse momenta) [19]. In the framework of the large-M t approximation, the NLO QCD corrections to the transversemomentum distribution of the SM Higgs boson were computed in Refs. [20][21][22][23]. Corrections to the large-M t approximation are considered in Ref. [24]. The numerical programs of Refs. [20,23] can also be used to evaluate arbitrary infrared-and collinear-safe observables up to NLO in the large-q T region and, in the case of Ref. [23], up to NNLO when q T = 0.
In the small-q T region (q T ≪ M H ), where the bulk of events is produced, the convergence of the fixed-order expansion is spoiled, since the coefficients of the perturbative series in α S (M 2 H ) are enhanced by powers of large logarithmic terms, ln m (M 2 H /q 2 T ). To obtain reliable perturbative predictions, these terms have to be resummed to all orders in α S . The method to systematically perform all-order resummation of classes of logarithmically-enhanced terms at small q T is known [25][26][27][28][29][30][31][32][33]. In the case of the SM Higgs boson, resummation has been explicitly worked out at leading logarithmic (LL), next-to-leading logarithmic (NLL) [34,35] and next-to-next-to-leading logarithmic (NNLL) [36] level.
The fixed-order and resummed approaches at small and large values of q T can then be matched at intermediate values of q T , to obtain QCD predictions for the entire range of transverse momenta. Phenomenological studies of the SM Higgs boson q T distribution have been performed in Refs. [37, 35, 38-42, 1, 43-46], by combining resummed and fixed-order perturbation theory at different levels of theoretical accuracy. A comparison of theoretical calculations [1,40,42,44] and of results from parton shower Monte Carlo generators [47][48][49][50] is presented in Ref. [51].
In the present paper we compute the Higgs boson q T distribution at the LHC by combining the most advanced perturbative information that is available at present: NNLL resummation at small q T and NLO perturbation theory at large q T . The first results of our calculation were presented in Refs. [1,52]. Here we perform a more complete phenomenological study and present a discussion of theoretical uncertainties.
The formalism used to obtain these results was briefly described in Refs. [33,1] and is illustrated in detail in the present paper. Three distinctive features are anticipated here. The resummation is performed at the level of the partonic cross section; this implies that the parton distributions are evaluated at the factorization scale µ F , which has to be chosen of the order of the hard scale M. The resummed terms are embodied in a form factor that is universal: it depends only on the flavour of the partons that initiate the hard-scattering subprocess at the Born level (e.g. qq annihilation in the case of Drell-Yan lepton pair production, and gg fusion in the case of Higgs boson production). A constraint of perturbative unitarity is imposed on the resummed terms, to the purpose of reducing the effect of unjustified higher-order contributions at large values of q T and, especially, at intermediate values of q T . The constraint implies that the total cross section at the nominal fixed-order accuracy (NLO or NNLO) is recovered upon integration over q T of the transverse-momentum spectrum (at NLL+LO or NNLL+NLO).
The paper is organized as follows. In Sect. 2 the resummation formalism is discussed in detail. After illustrating the general aspects of our approach in Sect. 2.1, we discuss the structure of the resummed cross section in Sect. 2.2. The relation to the standard b-space resummation is given in Sect. 2.3. Sect. 2.4 is devoted to the finite component of the cross section. In Sect. 3 we apply the resummation formalism to the production of the SM Higgs boson at the LHC. In Sect. 4 we draw our conclusions. In Appendix A we discuss the details of the exponentiation in the general multiflavour case. In Appendix B we illustrate the calculation of the Bessel integrals required in the computation of the perturbative expansion of the resummed cross section.

Transverse-momentum resummation
The formalism [1,33] that we use to compute the q T distribution of the Higgs boson applies to more general hard-scattering processes. Therefore, we describe it in general terms.

The resummation formalism: from small to large values of q T
We consider the inclusive hard-scattering process where the collision of the two hadrons h 1 and h 2 with momenta p 1 and p 2 produces the triggered final-state system F , accompanied by an arbitrary and undetected final state X. We denote by √ s the centre-of-mass energy of the colliding hadrons (s = (p 1 + p 2 ) 2 ≃ 2p 1 p 2 ). The observed final state F is a generic system of non-QCD partons such as one or more vector bosons (γ * , W, Z, . . . ), Higgs particles, Drell-Yan (DY) lepton pairs and so forth. We do not consider the production of strongly interacting particles (hadrons, jets, heavy quarks, ...), since in this case the resummation formalism of small-q T logarithms has not yet been fully developed.
Throughout the paper we limit ourselves to considering the case in which only the total invariant mass M and transverse momentum q T of the system F are measured. According to the QCD factorization theorem (see Ref. [53] and references therein), the corresponding transversemomentum differential cross section † dσ F /dq 2 T can be written as (2) † To be precise, when the system F is not a single on-shell particle of mass M , what we denote by dσ F /dq 2 T is actually the differential cross section M 2 dσ F /dM 2 dq 2 T .
where f a/h (x, µ 2 F ) (a = q f ,q f , g) are the parton densities of the colliding hadrons at the factorization scale µ F , dσ F ab /dq 2 T are the partonic cross sections,ŝ = x 1 x 2 s is the partonic centre-of-mass energy, and µ R is the renormalization scale. Throughout the paper 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 section is computable in QCD perturbation theory as a power series expansion in α S . We assume that at the parton level the system F is produced with vanishing q T (i.e. with no accompanying final-state radiation) in the lowest-order approximation, so that the corresponding cross section is dσ F cc /dq 2 T ∝ δ(q 2 T ). Since F is colourless, the lowest-order partonic subprocess, c +c → F , is either qq annihilation (c = q), as in the case of γ * , W and Z production, or gg fusion (c = g), as in the case of the production of the SM Higgs boson H.
As recalled in Sect. 1, higher-order perturbative contributions to the partonic cross section dσ F ab /dq 2 T contain logarithmic terms of the type ln m (M 2 /q 2 T ) that become large in the small-q T region (q T ≪ M). Therefore, we introduce the following decomposition of the partonic cross section in Eq. (2): The distinction between the two terms on the right-hand side is purely theoretical. The first term, dσ (res.) F ab , on the right-hand side contains all the logarithmically-enhanced contributions, (α n S /q 2 T ) ln m (M 2 /q 2 T ), at small q T , and has to be evaluated by resumming them to all orders in α S . The second term, dσ (fin.) F ab , is free of such contributions, and can be computed by fixed-order truncation of the perturbative series. More precisely, we define the 'finite' component dσ where the right-hand side vanishes order-by-order in perturbation theory. In particular, this implies that any perturbative contributions proportional to δ(q 2 T ) have been removed from dσ (fin.) F ab and included in dσ The 'resummed' component dσ F ab of the partonic cross section cannot, of course, be evaluated by computing all the logarithmic contributions in the perturbative series. However, as discussed in Sect. 2.2, these contributions can systematically be organized in classes of LL, NLL, . . . terms and, then, this logarithmic expansion can be truncated at a given logarithmic accuracy.
In summary, the q T distribution in Eq. (2) is evaluated, in practice, by replacing the partonic cross section on the right-hand side as follows The first and second terms on the right-hand side denote the truncation of the resummed and finite components at a given logarithmic accuracy and at a given fixed order, respectively. The resummed component gives the dominant contribution in the small-q T region, while the finite ‡ The notation X f.o. means that the quantity X is computed by truncating its perturbative expansion at a given fixed order in α S . starting from dσ ab f.o. , the usual perturbative series for the partonic cross section truncated at a given fixed order in α S , and subtracting from it the perturbative truncation of the resummed component at the same fixed order in α S : Moreover, we impose the condition: This matching procedure guarantees that the replacement in Eq. (5) retains the full information of the perturbative calculation up to the specified fixed order plus resummation of logarithmicallyenhanced contributions from higher orders. Equations (6) and (7) indeed imply that the matching is perturbatively exact, in the sense that the fixed-order truncation of the right-hand side of Eq. (5) exactly reproduces the customary fixed-order truncation of the partonic cross section in Eq.
(2). The (small-q T ) resummed and (large-q T ) fixed-order approaches are thus consistently combined without double-counting (or neglecting) of perturbative contributions and by avoiding the introduction of ad-hoc boundaries (such as, for instance, the choice of some intermediate value of q T as 'switching' point between the resummed and fixed-order calculations) between the large-q T and small-q T regions.
The resummed contributions that are present in the term dσ (res.) F ab l.a. of Eq. (5) are necessary and fully justified at small q T . Nonetheless they can lead to sizeable higher-order perturbative effects also at large q T , where the small-q T logarithmic approximation is not valid. To reduce the impact of these unjustified higher-order terms, we require that they give no contributions to the most basic quantity, namely the total cross section, that is not affected by small-q T logarithmic terms. We thus impose that the integral over q T of Eq. (5) exactly reproduces the fixed-order calculation of the total cross section. Since dσ (fin.) F ab is evaluated in fixed-order perturbation theory, the perturbative constraint on the total cross section is achieved by imposing the following condition: Equation (8) can be regarded, in some sense, as a unitarity constraint. As a matter of fact, the logarithmic contributions that are resummed in dσ (res.) F ab are, precisely speaking, plus distributions of the type [(α n S /q 2 T ) ln m (M 2 /q 2 T )] + . Therefore, it is quite natural to require that these resummed terms give a vanishing contribution to the total cross section. Note that the bulk of the q T distribution is in the region q T ∼ < M H . Since resummed and fixed-order perturbation theory controls the small-q T and large-q T regions respectively, the total cross section constraint mainly acts on the size of the higher-order contributions introduced in the intermediate-q T region by the matching procedure.
Another distinctive feature of the formalism illustrated so far is that we implement perturbative QCD resummation at the level of the partonic cross section. In the factorization formula (2), the parton densities are thus evaluated at the factorization scale µ F , as in the customary perturbative calculations at large q T . Although we are dealing with a process characterized by two distinct hard scales, q T and M, the dominant effects from the scale region q T ≪ M are explicitly taken into account through all-order resummation. Therefore, the central value of µ F and µ R has to be set equal to M H , the 'remaining' typical hard scale of the process. Then the theoretical accuracy of the resummed calculation can be investigated as in customary fixed-order calculations, by varying µ F and µ R around this central value.
At small values of q T , the perturbative QCD approach has to be supplemented with nonperturbative contributions, since they become relevant as q T decreases. A discussion on nonperturbative effects on the q T distribution of the SM Higgs boson is presented in Sect. 3.1.
The resummation and matching formalism, which we have so far illustrated in quite general terms, is set up to deal with the transverse-momentum region where q T ∼ < M. Resummation of small-q T logarithms cannot lead to any theoretical improvements in the large-q T region, where those logarithms are not the dominant contributions. When q T ∼ > M, the use of the resummation formalism is no longer justified (recommended), and we have to use the customary fixed-order perturbative expansion.

The resummed component
The method to systematically resum the logarithmically-enhanced contributions at small q T was set up [26][27][28][29][30] shortly after the first resummed calculation of the DY q T spectrum to double logarithmic accuracy [25]. The resummation procedure has to be carried out in the impactparameter space, to correctly take into account the kinematics constraint of transverse-momentum conservation. The resummed component of the transverse-momentum cross section in Eq. (3) is then obtained by performing the inverse Fourier (Bessel) transformation with respect to the impact parameter b. We write § dσ (res.) F ab where J 0 (x) is the 0th-order Bessel function.
The perturbative and process-dependent factor W F ab embodies the all-order dependence on the large logarithms ln M 2 b 2 at large b, which correspond to the q T -space terms ln M 2 /q 2 T that are logarithmically enhanced at small q T (the limit q T ≪ M corresponds to Mb ≫ 1, since b is the variable conjugate to q T ). Resummation of these large logarithms is better expressed by defining the N-moments ¶ W N of W with respect to z = M 2 /ŝ at fixed M: The resummation structure of W F ab, N can indeed be organized in exponential form, as discussed below.
In the following of this subsection, the subscripts denoting the flavour indices are understood. More precisely, we present the resummation formulae in a simplified form, which is valid when there is a single species of partons. This simplified form illustrates more clearly the key structure and the main features of the resummed partonic cross section. The generalization to considering more species of partons does not require any further conceptual steps: it just involves algebraic complications, which are discussed in Sect. 2.3 and in Appendix A.
The logarithmic terms embodied in W F ab, N are due to final-state radiation of partons that are soft and/or collinear to the incoming partons. Their all-order resummation can be organized [33] in close analogy to the cases of soft-gluon resummed calculations for hadronic event shapes in hard-scattering processes [54][55][56][57] and for threshold contributions to hadronic cross sections [58,59]. We write The function H F N does not depend on the impact parameter b and, therefore, it contains all the perturbative terms that behave as constants in the limit b → ∞. The function G includes the complete dependence on b and, in particular, it contains all the terms that order-by-order in α S are logarithmically divergent when b → ∞. This factorization between constant and logarithmic terms involves some degree of arbitrariness [56], since the argument of the large logarithms can always be rescaled as ln M 2 b 2 = ln Q 2 b 2 + ln M 2 /Q 2 , provided that Q is independent of b and that ln M 2 /Q 2 = O(1) when bM ≫ 1. To parametrize this arbitrariness, on the right-hand side of Eq. (12) we have introduced the scale Q, such that Q ∼ M, and we have defined the large logarithmic expansion parameter, L, as where the coefficient b 0 = 2e −γ E (γ E = 0.5772 . . . is the Euler number) has a kinematical origin (the use of b 0 in Eq. (13) in purely conventional: it simplifies the algebraic expression of G).
The role played by the auxiliary scale Q (which we name the 'resummation scale') in the context of the resummation program is analogous to the role played by the renormalization (factorization) scale in the context of renormalization (factorization). Although the resummed cross section W F N does not depend on Q when evaluated to all perturbative orders, its explicit dependence on Q appears when W F N is computed by truncation of the resummed expression at some level of logarithmic accuracy (see below). As in the case of µ R and µ F , we should set Q at the central value Q = M; variations of the resummation scale Q around this central value can then be used to estimate the uncertainty from yet uncalculated logarithmic corrections at higher orders. Note that the resummation scale dependence of W F N should not be confused with the 'resummation scheme' dependence considered in Ref. [33]. In fact, as shown in Sect. 2.3, W F N is exactly independent of the resummation scheme.
All the large logarithmic terms α n S L m with 1 ≤ m ≤ 2n are included in the form factor exp{G}. More importantly, all the logarithmic contributions to G with n + 2 ≤ m ≤ 2n are vanishing. This property, which is called exponentiation, follows [26][27][28][29][30] from the perturbative dynamics of (abelian and non-abelian) gauge theories and from kinematics factorization in impact parameter space. Thus, the exponent G can systematically be expanded as where α S = α S (µ 2 R ) and the functions g (n) (α S L) are defined such that g (n) = 0 when α S L = 0. Thus the term Lg (1) collects the LL contributions α n S L n+1 ; the function g (2) resums the NLL contributions α n S L n ; g (3) controls the NNLL terms α n S L n−1 , and so forth. Note that in the context of the resummation approach, the parameter α S L is formally considered as being of order unity. Thus, the ratio of two successive terms in the expansion (14) is formally of O(α S ) (with no L enhancement). In this respect, the resummed logarithmic expansion in Eq. (14) is as systematic as any customary fixed-order expansion in powers of α S .
The function H F N in Eq. (12) does not contain large logarithmic terms to be resummed. It can be expanded in powers of α S = α S (µ 2 R ) as is the lowest-order partonic cross section for the hard-scattering process in Eq. (1).
Two other general aspects of the resummed partonic cross section W F N are the factorization scale (and scheme) dependence and the process dependence. As discussed below, the form factor exp{G} does not depend on both the factorization scale (and scheme) and the specific hardscattering process.
The hadronic cross section on the left-hand side of Eq. (2) is a physical observable and cannot depend on the factorization scale µ F . In practice, the evaluation of the right-hand side at a certain perturbative accuracy introduces the µ F dependence of the partonic cross section dσ F ab . This dependence is perturbatively balanced by the µ F dependence of the parton densities f a/h (x, µ 2 F ). Note that the parton densities in Eq. (2) do not depend on the transverse momentum q T (or on the impact parameter b). Recall also that we implement transverse-momentum resummation at the level of the partonic cross section dσ F ab , by using Eqs. (9) and (12). Therefore, any µ F dependence of the parton densities cannot introduce any logarithmic dependence on b in the form factor exp{G}. In other words, the perturbative expansion (15) of the function H F N depends on µ F , while the exponent G of the form factor and its corresponding logarithmic functions g (n) N in Eq. (14) do not depend on µ F and on the factorization scheme used to define the parton densities.
As explicitly shown in Sect. 2.3, the form factor exp{G} in Eq. (12) does not depend on the final-state system F produced in the hard-scattering process of Eq. (1). The form factor is process independent: it is produced by universal soft and collinear radiation from the QCD partons entering the hard-scattering process (when the simplification of considering a single parton species is removed, there are various process-independent form factors for the various partonic channels). The dependence on the process is fully taken into account by the hard-scattering function H F N , which embodies contributions produced by virtual corrections at transverse-momentum scales q T ∼ M.
The truncation W F N l.a. of the resummed cross section at a given logarithmic accuracy is defined as follows. At LL accuracy, we include the function g (1) in the exponent G and we approximate H F N by the Born cross section σ (0) F . At NLL accuracy, we include the functions g (1) and g (2) N and the coefficient H . At NNLL accuracy, we also include g N at NLL accuracy is that the combined effect of α S H F (1) N and Lg (1) (α S L) leads to logarithmic contributions, α n S L n , that are of the same order as those in g  The logarithmic truncation of the resummed component of the cross section can then be combined, as in Eq. (5), with the fixed-order expansion of the finite component in Eq. (6). The NLL+LO result is obtained by supplementing NLL resummation with the LO expansion at large q T . The NNLL+NLO result combines NNLL resummation with the NLO expansion at large q T . This procedure for combining the resummed and fixed-order approaches exactly satisfies the matching conditions in Eqs. (4) and (7). Note that the fulfilment of the matching conditions is not completely trivial. For instance, if H To reduce the impact of unjustified resummed logarithms in the large-q T region, we use a procedure inspired by that introduced in Ref. [55] to deal with kinematical constraints when performing soft-gluon resummation in e + e − event shapes. We consider the exponent G(α S , L) of the form factor in Eqs. (12) and (14) and we perform the replacement In other words, in the argument of G(α S , L) we replace the logarithmic variable L with the variable L defined as Comparing the definitions in Eqs. (13) and (17), we see that in the resummation region Qb ≫ 1 we have L = L + O(1/(Qb) 2 ), and thus the replacement in Eq. (16) is fully legitimate * * to arbitrary logarithmic accuracy. Although the variables L and L are equivalent to organize the We recall that there is a mismatch of notation between the q T distribution at q T ∼ M and the total cross section. The LO (NLO) term of the finite component of the q T distribution contributes to the total cross section at NLO (NNLO). * * Note that the replacement in Eq. (16) introduces an explicit dependence of dσ (res.) F on the resummation scale Q. Owing to the matching procedure in Eq. (6), this dependence is balanced by the Q dependence of the dσ (fin.) F . resummation formalism in the region Qb ≫ 1, they lead to a different behaviour of the form factor at small values of b (i.e. large values of q T ): when Qb ≪ 1, we have L → 0 and exp{G(α S , L)} → 1. Therefore, performing the replacement in Eq. (16), we reduce the effect produced by the resummed contributions in the small-b region, where the use of the large-b resummation approach is not justified.
In particular, since exp{G(α S , L)} = 1 at b = 0, using Eqs. (9) and (12) we obtain the relation which simply follows from the fact that the value at b = 0 of the (b-space) Fourier transformation of the q T distribution is equal to the integral over q T of the q T distribution itself. Since the hard cross section H F is evaluated in fixed-order perturbation theory, the relation (18) implies that the replacement in Eq. (16) also allows us to implement the perturbative constraint (8) on the total cross section. More precisely, the integral over q T of the q T distribution dσ F /dq T at NLL+LO (NNLL+NLO) accuracy exactly reproduces the calculation of the total cross section at NLO (NNLO).
The purpose of the transverse-momentum resummation program [26][27][28][29][30] is to explicitly evaluate the logarithmic functions g (n) N of Eq. (14) in terms of few coefficients that are perturbatively computable. As illustrated in Sect. 2.3, this goal is achieved by showing that the all-order resummation formula (14) has the following integral representation: where A(α S ) and B N (α S ) are perturbative functions The coefficients A (n) and B   21)), which can be extracted from the logarithmic terms in the perturbative expansion of the q T distribution at the n-th relative order in α S . Therefore, the customary fixed-order computation of the q T distribution is sufficient to obtain the full information that is necessary to explicitly perform resummation at the required logarithmic accuracy.
By inspection of the q 2 integration in Eq. (19), it is evident that the exponent G N of the processindependent form factor in Eq. (12) has the logarithmic structure of Eq. (14). The functions g (n) N depend on the coefficients in Eqs. (20) and (21), and the functional dependence is completely specified by Eq. (19). More precisely (see Eqs. (22)-(24)), the LL function g N depends also on A (2) and B (1) N , the NNLL function g (3) N depends also on A (3) and B (2) N , and so forth. Starting from the integral representation in Eq. (19), the explicit functional form of the functions g N (for arbitrary values of n) can easily be computed by using the method that is described in Appendix C of Ref. [17].
The LL, NLL and NNLL functions g N have the following explicit expressions † † : where and β n are the coefficients of the QCD β function: (27) † † Note that the functional form of the functions g (n) N is exactly the same as that of the functions that appear in the calculation of the energy-energy correlation in e + e − annihilation [60].
The explicit expression of the first three coefficients, β 0 , β 1 and β 2 , is [61] where N f is the number of QCD massless flavours and the SU(N c ) colour factors are Note that the functions g  This type of singularities * is a common feature of all-order resummation formulae of soft-gluon contributions. Within a perturbative framework, these singularities have to be regularized. A possible regularization procedure consists in introducing a 'minimal prescription', such as those introduced in Refs. [59] (in the case of threshold resummation) and [62,44] (in the case of b-space or joint resummation). In the case of b-space resummation, other procedures are to use the 'b * prescription' of Ref. [29], by freezing the integration over b below a fixed upper limit, or more simply, to introduce a cut-off at a very large (but smaller than b L ) value of b [63]. Admittedly, when the non-perturbative contributions are sizeable, they have to be properly included, according to the prescription used to regularize the singularities.

Sudakov form factor, universal form factor and perturbative coefficients
The b-space resummation approach was fully formalized by Collins, Soper and Sterman [28,32] in terms of perturbative coefficients. Considering the generic hard-scattering process in Eq. (1), the transverse-momentum differential cross section in Eq.
(2) is written as where the dots on the right-hand side stand for terms that are not logarithmically enhanced at small q T (large b). Note that Eq. (29) regards the hadronic cross section (and not the partonic cross section in Eq. (10)). Therefore, the b-space function W F (b, M, s), which embodies the logarithmically-enhanced terms, depends on the parton densities of the colliding hadrons. The all-order resummation of the large logarithms ln(M 2 b 2 ) in the region Mb ≫ 1 is accomplished by showing that the N-moments W N (b, M) of W (b, M, s) with respect to z = M 2 /s at fixed M can be recast in the following form [32,33]: where f a/h, N (µ 2 ) are the N-moments of the parton density f a/h (z, µ 2 ), and σ (0) cc, F is the lowestorder cross section for the partonic subprocess c +c → F . The function S c (M, b) is the Sudakov form factor of the quark (c = q,q) or of the gluon (c = g), and it has the following expression † : .
The functions A, B, C and H F in Eqs. (30) and (31) are perturbative series in α S : The functions A c , B c and C ab are process independent, while H F c depends on the specific hardscattering process.
The resummation formulae (30) and (31) are invariant under the following 'resummation scheme' transformations [33]: The invariance can easily be proven by using the following renormalization-group identity (see Eq. (27)): which is valid for any perturbative function h(α S ). † In Ref. [32] the upper limit of the integral in Eq. (31) is set to C 2 M 2 , where C 2 is an arbitrary factor. The scale C 2 M 2 is thus related to the resummation scale Q 2 in Eq. (19).
The physical origin of the resummation scheme invariance of Eq. (30) is discussed in Ref. [33]. The invariance implies that the factors H F c , S c (more precisely, the function B c ) and C ab are not unambiguously computable order by order in perturbation theory. In other words, these factors can be unambiguously defined only by choosing a 'resummation scheme'. The choice of a resummation scheme amounts to defining H F c (or C ab ) for a single process. More precisely, H F c has to be defined for two processes: one process that is controlled, at the lowest perturbative order, by qq annihilation (c = q,q) and another process that is controlled by gg fusion (c = g). Having done that, the process-dependent factor H F c and the universal (process-independent) factors S c and C ab are unambiguously determined for any other process of the type in Eq. (1).
Such a form is certainly consistent since, by choosing h(α S ) = H F c (α S ) and using the invariance under the transformation in Eq. (36), it is always possible to set H F c (α S ) = 1 on a process-dependent basis. Note that this procedure does not correspond to the definition of a resummation scheme. Indeed, the corresponding Sudakov form factor S F c and the functions C F ab turn out to be process-dependent quantities, as pointed out by the explicit and general calculation of B (2) c and C (1) ab (z) in Ref. [36]. For example, in the case of gg fusion processes, the Sudakov form factors for the production of a scalar and a pseudoscalar Higgs boson turn out to be different and to have even a different dependence on the mass of the top quark.
Comparing the partonic and the hadronic cross sections in Eqs. (10) and (30), we see that the resummed factors W F ab and To express the resummed partonic cross section W F ab in terms of the perturbative coefficients in Eqs. (32)-(35), we have to use Eq. (30) and substitute the parton densities f a/h, N (b 2 0 /b 2 ) for the same parton densities evaluated at the factorization scale µ F . The substitution can be done by using where the QCD evolution operator U ab, N (µ 2 , µ 2 0 ) fulfils the evolution equations and γ ab, N (α S ) are the parton anomalous dimensions or, more precisely, the N-moments of the customary Altarelli-Parisi splitting functions P ab (α S , z) [64]: We finally obtain [33] which relates the resummed partonic cross section in Eq. (10) to the perturbative coefficients in Eqs. (32)- (35) and the anomalous dimensions coefficients in Eq. (41).
In the following we explicitly show how Eq. (42) is related to the exponential structure of Eq. (12) in the case with a single species of partons. The general case with partons of different flavours is discussed in Appendix A. Here we only anticipate that the generalization of Eq. (12) to the multiflavour case ‡ simply involves a sum of exponential terms, namely where the index {I} labels a set of flavour indices (which is precisely specified in Appendix A).
Within the simplified treatment in which there is a single species of partons, the resummed partonic cross section in Eq. (42) can easily be recast in the factorized exponential form of Eqs. (12) and (19). To this aim, we first use the identity (37) ). Then, we insert in Eq. (42) the solution of the evolution equation (40): We finally obtain the exponential form in Eq. (19), where the perturbative function A(α S ) is exactly the perturbative function in Eq. (32), and the function B N (α S ) is given as follows in terms of the perturbative functions in Eqs. (27), (33), (34) and (41): The expression of the hard-process function H F N in Eq. (12) is Note that, as discussed in Sect. 2.2, the form factor exp{G} and, hence, the perturbative functions A(α S ) and B N (α S ) in Eq. (19) do not depend on the factorization scale µ F . As a consequence, the functions A(α S ) and B N (α S ) are also independent of the factorization scheme used to define the parton densities. Since, as is well known, the anomalous dimensions γ ab, N (α S ) do depend on the factorization scheme, the relation (45) implies that both the perturbative functions B c (α S ) and C ab (α S ) depend on the factorization scheme in such a way that B N (α S ) turns out to be factorization-scheme independent.
As anticipated in Sect. 2.2, the form factor exp{G} does not depend on the final-state system F produced in the hard-scattering process. From Eqs. (19) and (45), this independence is a simple consequence of the process independence of each of the perturbative functions A c (α S ), B c (α S ), C ab (α S ) and γ ab, N (α S ).
The relation (45) also implies that the form factor exp{G} does not depend on the resummation scheme used to express the various factors in the resummation formulae (30) and (31) (we recall that the customary Sudakov form factor S c (M, b) in Eq. (31) does instead depend on the resummation scheme). It is indeed straightforward to show that the function B N (α S ) in Eq. (45) is invariant under the resummation-scheme transformations in Eq. (36).
Unlike the form factor exp{G}, the non-logarithmic function H F N in Eq. (46) explicitly depends on the factorization scale µ F , on the factorization scheme (through C ab, N (α S ) and γ ab, N (α S )) and on the final-state system F (through σ The universal (i.e. independent of the process and of the factorization and resummation schemes) perturbative function A c (α S ) in Eqs. (20) and (32) where is not yet known. In our quantitative study of transverse-momentum resummation at NNLL accuracy (see Sect. 3), we assume that the value of A (2) c in Eq. (47) are exactly equal to those of the related perturbative function that controls threshold resummation [58] in the MS factorization scheme. Note, however, that the two soft-gluon functions A c (α S ) do not necessarily coincide at high perturbative orders since, for instance, the soft-gluon function for transverse-momentum resummation is universal while the soft-gluon function for threshold resummation depends on the factorization scheme.
The first-order coefficient B (1) c, N of the universal perturbative function B N (α S ) in Eqs. (21) and (45) is B (1) where [30,34] Note that, since the LO anomalous dimensions γ (1) cc, N are universal, the NLL coefficients B cc, N + 2γ (2) or, equivalently, by performing the inverse Mellin transformation to z-space: The value of the quark coefficient B (2) q can be obtained by using the results of Ref. [67] for the coefficients B (2) q and C (1) qq (z) of the DY process. These results are confirmed by the general (process-independent) calculation of Ref. [36], which considers both the qq-annihilation and the gluon fusion channels. From the results of Ref. [36] we obtain the value of the gluon coefficient B (2) g , and we can also explicitly check the universality of both B (2) q and B (2) g . To write down the expression of B (2) c , we recall that the second-order term P (2) cc (z) of the Altarelli-Parisi splitting functions P cc (α S , z) has the following general dependence on z: where A (2) c is the coefficient in Eq. (47), 1/(1 − z) + is the customary 'plus'-distribution and P (2)reg cc (z) denotes all the remaining and less singular (when z → 1) contributions to P The first-order coefficients C qg and C gq in Eq. (34) do not depend on the process and on the resummation scheme, and were first computed in Refs. [67] and [35], respectively. Their expressions in the MS factorization scheme are The flavour-diagonal first-order coefficients C where we have defined The coefficient H depends on the process and is explicitly known for several processes (see Ref. [36] and references therein).
To complete the resummation program at NNLL, the coefficient H is also needed. This coefficient is not known in analytic form for any hard-scattering process. Nonetheless, within our resummation formalism, it can be determined for any hard-scattering process whose corresponding total cross section is known at NNLO. This point is discussed in detail at the end of Sect. 2.4. F ab /dq 2 T does not contain any perturbative contributions proportional to δ(q 2 T ) (these contributions and all the logarithmicallyenhanced terms at small q T are included in dσ

The finite component
where the subscript LO (NLO) denotes the perturbative truncation of the various cross sections at the leading order (next-to-leading order) in the region where q T = 0. The extension of Eqs. (58) and (59) at still higher perturbative order is straightforward.
The contributions dσ F ab f.o. on the right-hand side of Eqs. (58) and (59)  in Eqs. (10). To this purpose, we define the perturbative coefficients Σ (n) as follows: cc, F (M) and, in general, the power p cF depends on the lowest-order partonic subprocess c +c → F . In Eq. (60), W F ab is the resummed cross section on the right-hand side of Eq. (10). Note, however, that Eq. (60) depends on the resummation scale Q 2 . The dependence on the resummation scale has been introduced in Eqs. (10) and (12)  The perturbative expansion of Eq. (12) or, more precisely, of Eq. (42) gives where the dependence on the scale ratios M 2 /µ 2 R , M 2 /µ 2 F and M 2 /Q 2 is understood. The extension of Eqs. (61) and (62) to the higher order terms Σ F (n) cc←ab (z, L) with n ≥ 3, is straightforward. The b-independent coefficients Σ F (1;k) (z), H F (1) (z), Σ F (2;k) (z) and H F (2) (z) are more easily presented by considering their N-moments with respect to the variable z. We have The right-hand side of Eqs. (63)-(70) is expressed in terms of the resummation-scheme independent coefficients given in Sect. 2.3 and of the logarithms ℓ Q , ℓ F and ℓ R defined in Eq. (57). To explicitly exhibit the independence of the resummation scheme we can, for example, rewrite the contribution in the third line of Eq. (69) in terms of the resummation-scheme independent coefficients B (2) c N (see Eq. (50)) and C (1) ab, N with a = b (see Eq. (55)): Inserting Eqs. (60)- (62) in Eq. (10), performing the integral over the impact parameter b, and removing the contributions proportional to δ(q 2 T ) (for example, all the contributions coming from H F (n) cc←ab in Eq. (60)), we obtain the following expressions for the fixed-order contributions dσ (res.) F ab f.o. on the right-hand side of Eqs. (58) and (59): On the right-hand side of Eqs. (72) and (73), the dependence on q T is fully embodied in the functions I n (q T /Q), which are obtained by the following Bessel transformation: The term ln n (1 + Q 2 b 2 /b 2 0 ) = L n in the integrand comes from the replacement L → L (see Eq. (16)). In customary implementations of b-space resummation, one has to consider the Bessel transformation of powers of ln n (Q 2 b 2 /b 2 0 ) = L n , which can be expressed in terms of powers of ln n (Q 2 /q 2 T ). The functions I n (q T /Q) have instead a more involved functional dependence on q T . As shown in Appendix B, this functional dependence can be expressed in terms of K ν (q T /Q), the modified Bessel function of imaginary argument that is defined by the following integral representation: We conclude this section with some observations on the hard-scattering function H F cc←ab . This function is resummation-scheme independent, but it depends on the specific hard-scattering subprocess c +c → F . The coefficients H F (n) cc←ab of its perturbative expansion can be determined by performing a customary perturbative calculation of the q T distribution in the limit q T → 0. Moreover, as discussed in Sect. 2.2, within our resummation formalism H F controls the strict perturbative normalization of the corresponding total cross section (i.e. the integral of the q T distribution). This property can be exploited to determine the coefficients H F (n) cc←ab in a different manner, that is, from the perturbative calculation of the total cross section.
To illustrate this point we consider the total cross section,σ tot F ab , at the partonic level, and we evaluate the q T spectrum on right-hand side according to the decomposition in terms of 'resummed' and 'finite' components (see Eq. (3)). Then we use Eq. (18) to integrate the resummed component over q T , and we obtain This expression is valid order by order in QCD perturbation theory. Once the perturbative coefficients of the fixed-order expansions ofσ tot F ab , H F ab and dσ (fin.) F ab /dq 2 T are all known, the relation (77) has to be regarded as an identity, which can explicitly be checked. Note, however, that since the fixed-order truncation dσ F ab /dq 2 T at N n−1 LO, without the need of explicitly computing the small-q T behaviour of the spectrum dσ F ab /dq 2 T at N n LO. For example, at NLO Eq. (77) gives where α S = α S (µ 2 R ) and we have used At NNLO Eq. (77) gives and the generalization at still higher orders n (n > 2) is In our study of the transverse-momentum spectrum of the Higgs boson at NNLL accuracy (see Sect. 3), we use Eq. (80) to obtain a numerical value for the corresponding perturbative coefficient H (2) .

The q T spectrum of the Higgs boson at the LHC
In this section we apply the resummation formalism described in Sect. 2 to the production of the SM Higgs boson at the LHC.
We consider the gluon fusion production mechanism gg → H, whose Born level cross section in Eqs. (15) and (60) is As recalled in Sect. 1, this implementation of the large-M t approximation is expected to produce an uncertainty that is smaller than the uncertainties from yet uncalculated perturbative terms from higher orders.
We compute the Higgs boson differential cross section dσ/dq T at the LHC (pp collisions at √ s = 14 TeV) and present quantitative results at NLL+LO and NNLL+NLO accuracy.
As discussed in Sect. 2.2, at NLL+LO accuracy the resummed component in Eq. (12) is evaluated by including the functions g (1) and g (2) N in Eq. (14) and the coefficient H in Eq. (15), and then it is matched with the fixed-order contribution evaluated at the LO (i.e. at O(α 3 S )) in the large-q T region. The functions g (1) and g (2) N are process independent and given in terms of the universal coefficients A (1) , A (2) and B where the coefficient C (1) gq,N is the Mellin transformation of Eq. (55). Of course, these processindependent coefficients are exact, i.e. they are not affected by the large-M t approximation. The flavour diagonal coefficient H gg←gg,N is instead process dependent; therefore it depends on M t and, in the large-M t approximation, it is given by [39,35] H H(1) where, for simplicity, the scale-dependent terms have been dropped (i.e. we have set µ R = µ F = Q = M H in Eq. (65)).
At NNLL+NLO accuracy the function g  70)) is not known in analytic form. We thus exploit Eq. (80), which follows from the constraint of perturbative unitarity, to extract the numerical value of H H(2) N from the knowledge of the total cross section at the NNLO [14]. The scale-independent part of H F (2) gg←gg, N can be written as where the M t -dependent contribution on the right-hand side is obtained from the results in Refs. [16,68], and c N does not depend on M t in the large-M t approximation. Since from Eq. (84) we know that C (1) gg, N is actually independent of N, the N dependence of c N can only follow from that of C cc←ab, N can numerically be neglected, and that the coefficient c N in Eq. (86) can numerically be approximated by an N-independent value, c N ≃ 178. 75. This numerical approximation is pretty good, since the integral of the NNLL+NLO spectrum reproduces the NNLO total cross section to better than 1% accuracy in a wide Higgs mass range, 100 GeV ∼ < M H ∼ < 300 GeV, at the LHC.
We recall that the functions g  24)). The singular behaviour is related to the presence of the Landau pole in the perturbative running of the QCD coupling α S (q 2 ). As mentioned at the end of Sect. 2.2, a practical implementation of the resummation procedure requires a prescription to deal with these singularities. In our numerical study we follow Ref. [62] and deform the integration contour in the complex b space. In particular we choose the two integration branches as We have checked that the result is very mildly dependent on the choice of φ. We have also used the simpler procedure of integrating over the real b-axis, using a sharp cut-off at a large value of b, and checking the independence of the actual value of the cut-off. We found that the numerical differences between the results obtained by these two procedures are negligible.
Our complete calculation of the q T spectrum of the Higgs boson at the LHC is implemented in the numerical code HqT, which can be downloaded from [69] together with some accompanying notes. This code is a slightly modified and numerically improved version of the code used in Ref. [1]: the most important difference regards the computation of the finite component. In Ref. [1] we used the Monte Carlo program of Ref. [20] to compute the fixed-order contribution to the q T cross section at LO and NLO. Here we have implemented the analytic calculation of Glosser and Schmidt [22]. Although the two methods are in principle equivalent, the use of the analytic calculation allows us to achieve a faster numerical stability in the small-q T region. In the next subsection we present a selection of numerical results that can be obtained with our code. We also include a discussion of theoretical uncertainties.

Numerical results at the LHC
To compute the hadronic cross section, we use the MRST2004 set [70] of parton distribution functions. As for the perturbative order of the parton densities and α S , at variance with Ref. [1], we adopt here the following choice. At NLL+LO we use NLO parton densities and 2-loop α S , whereas at NNLL+NLO we use NNLO parton densities and 3-loop α S . This choice is perfectly consistent in the small q T region, since the corresponding partonic cross section is dominated by the resummed component evaluated at NLL and NNLL accuracy, respectively. The choice is fully justified also at intermediate values of q T , where the calculation of the partonic cross section is driven by the small-q T resummation and strongly constrained by the total cross section at NLO and NNLO, respectively. At large values of q T , q T ∼ M H , our evaluation of the partonic cross section is dominated by the fixed-order contributions at LO and NLO, respectively. Therefore, our choice introduces a formal mismatch with respect to the customary use of parton densities and α S . However, as shown and discussed later in this subsection, this formal mismatch does not lead to any inconsistencies at the quantitative level.
The NLL+LO spectrum with M H = 125 GeV is shown in Fig. 1. In the left-hand side, the full NLL+LO result (solid line) is compared with the LO one (dashed line) at the default scales µ F = µ R = Q = M H . We see that the LO calculation diverges to +∞ as q T → 0. The effect of the resummation, which is relevant below q T ∼ 100 GeV, leads to a physically well-behaved distribution: it has a kinematical peak at q T ∼ 12 GeV and vanishes as q T → 0. The LO finite component of the spectrum (dotted line), which is defined in Eq. (58), is also shown: as expected it dominates when q T ∼ M H and vanishes as q T → 0. Note, however, that the contribution of the finite component is sizeable in the intermediate-q T region (about 20% at q T ∼ 50 GeV) and not yet negligible at small values of q T (about 8% around the peak region). This underlies the importance of a careful and consistent matching between the resummed and fixed-order calculations. In the right-hand side of Fig. 1 we show the NLL+LO band as obtained by varying µ F and µ R simultaneously and independently in the range 0.5M H ≤ µ F , µ R ≤ 2M H with the constraint 0.5 ≤ µ F /µ R ≤ 2 (the resummation scale is kept fixed at Q = M H ). The scale dependence increases from about ±15% at the peak to about ±20% at q T = 100 GeV. The integral over q T of the NLL+LO spectrum is in agreement with the value of the NLO total cross section to better than 1%, thus proving the numerical accuracy of the code.
The NNLL+NLO results at the LHC are shown in Fig. 2. In the left-hand side, the full result (solid line) is compared with the NLO one (dashed line) at the default scales µ F = µ R = Q = M H . The NLO result diverges to −∞ as q T → 0 and, at small values of q T , it has an unphysical peak (the top of the peak is above the vertical scale of the plot) that is produced by the numerical compensation of negative leading logarithmic and positive subleading logarithmic contributions. The resummed result is physically well-behaved at small q T . The NLO finite component of the spectrum (dotted line), which is defined in Eq. (59), vanishes smoothly as q T → 0; its contribution amounts to about 10% in the peak region, about 17% at q T ∼ 25 GeV and about 35% at q T ∼ 50 GeV. This shows both the quality and the relevance of the matching procedure.
We find that the contribution of A (3) (recall from Sect. 2.3 that we are using an educated guess The right-hand side of Fig. 2 shows the scale dependence computed as in Fig. 1. The scale dependence is now about 8% at the peak and increases to about 20% at q T = 100 GeV. To better illustrate the main features of the dependence on the scales µ R and µ F , we present numerical results at two fixed values of q T in Figs. 3 and 4. In Fig. 3  As expected from the QCD running of α S , the cross sections typically decrease when µ R increases around the characteristic hard scale M H , at fixed µ F = M H . In the case of variations of µ F at fixed µ R = M H , we observe the opposite behaviour. This is not unexpected, since when M H = 125 GeV the cross section is mainly sensitive to partons with momentum fraction x ∼ 10 −2 , and in this x-range scaling violations of the parton densities are (moderately) positive. Varying the two scales simultaneously (µ F = µ R ) leads to a partial compensation of the two different behaviours. As a result, the scale dependence is mostly driven by the renormalization scale, because the lowest-order contribution to the process is proportional to α 3 S , a (relatively) high power of α S .
Comparing the LO with the NLL+LO results and the NLO with the NNLL+NLO results, we see that the scale dependence of the resummed results (solid lines) is smaller than that of the corresponding fixed-order results (dashed lines): the LO and NLL+LO curves have a comparable slope, but the NLL+LO results are higher; the NLO and NNLL+NLO results have smaller differences, but the slope of the NNLL+NLO curve is flatter. In summary, resummation reduces the scale dependence of the fixed-order calculations also in the region of intermediate values of q T .
In Fig. 4 we report analogous results at a smaller value of q T , namely q T = 15 GeV. The qualitative behaviour is similar to the one in Fig. 3. In this region of small transverse momenta the fixed-order result is no longer reliable (see Figs. 1 and 2), but its relative scale dependence does not increase and is even smaller than at q T = 50 GeV. This is due to the fact that the fixed-order cross section is much larger than at higher values of q T . The slope of the resummed results (solid lines) is sizeably flatter than that of the corresponding fixed-order results (dashed lines). We also notice a slight reduction in the scale dependence of the resummed results compared to Fig. 3, especially at NNLL+NLO accuracy.
In Fig. 5 the NLL+LO and NNLL+NLO bands shown in Figs. 1 and 2 are compared. We see that the NNLL+NLO band (solid lines) is smaller than the NLL+LO one (dashed lines) and overlaps with the latter at q T ∼ < 100 GeV. This suggests a good convergence of the resummed perturbative expansion. This result is confirmed by the inset plot, that shows the NNLL+NLO band normalized to the NLL+LO result at central value of the scales. This q T -dependent K factor, is stable, around the values 1.1-1.2, in the central region of the inset plot, and it increases (de-creases) drastically when q T ∼ > 50 GeV (q T ∼ < 2 GeV). In the large-q T region, the effect of perturbative higher-order corrections is known to be important [20][21][22]. At very small values of q T , non-perturbative effects are definitely expected to be relevant. We observe that a naive rescaling of the NLL+LO result by a constant (i.e. independent of q T ) K factor would not reproduce the NNLL+NLO result over the entire q T -range. The nice convergence of the resummed perturbative expansion suggested by Fig. 5 should be contrasted with the results in Fig. 6, where the corresponding fixed-order bands, computed as in Fig. 5, are shown. The results in Fig. 6 have no physical significance in the small-q T region, owing to the non-convergence of the fixed-order expansion herein. When q T ∼ > 25 GeV, we see that the scale dependence of the NLO (LO) result is larger than the one of the corresponding NNLL+NLO (NLL+LO) result in Fig. 5. More importantly, we see that the LO and NLO bands do not overlap. This implies that the scale dependence enclosed by these bands certainly underestimates the true theoretical uncertainty from missing higher-order terms. Equivalently, we can say that the uncertainty of these fixed-order calculations is more reliably estimated by performing scale variations over a range of scales that is wider than that used in Fig. 6. All this indicates a poor convergence of the fixed-order perturbative expansion at intermediate values of q T As mentioned at the beginning of this subsection, in our resummed calculations at NLL+LO and NNLL+NLO accuracy we use parton densities and α S at perturbative orders that are different from those customarily used in fixed-order calculations at LO and NLO, respectively. Indeed, the consistent procedure at large values of q T would be to use LO densities with 1-loop α S at the LO, and NLO densities with 2-loop α S at the NLO. We have also explained why our procedure is justified in the intermediate-q T region, and we have postponed the discussion on the large-q T region. To come back to this point, in Fig. 7 we compare our NLL+LO and NNLL+NLO results with the customary NLO results, which are obtained by using NLO parton densities and 2-loop α S . We also include the corresponding bands, computed from scale variations. In the left-hand side we see that in the intermediate-q T region our NLL+LO result catches the bulk of the NLO effect. Obviously, at large q T , the inclusion of NLO corrections is necessary. In the right-hand side, the calculations at NNLL+NLO accuracy and at the NLO are compared. In spite of the fact that the two calculations use different parton densities and α S , the corresponding bands show a very good overlap when q T ∼ M H . We thus conclude that, within the NLO theoretical uncertainty, the two calculations are perfectly compatible at the quantitative level in the large-q T region, q T ∼ M H .
In Fig. 8 (Fig. 9)  We first comment on the behaviour at large transverse momenta, which is best visible looking at the plots on the right of Figs. 8 and 9. We see that the NLL+LO cross section can become negative if Q = 2M H . This behaviour should not be regarded as particularly worrisome: it takes place when q T > M H , where the use of the resummation formalism is not anymore justified. In general, the cross section has a better behaviour at large q T when the resummation scale has the values Q = M H , M H /2, M H /4. In particular, at large-q T the results of the fixed-order calculation at LO (NLO) accuracy are very well approximated by the NLL+LO (NNLL+NLO) calculation with Q = M H /2; the line corresponding to the LO (NLO) results is not shown in the plot on the right of Fig. 8 (Fig. 9), since it is hardly distinguishable from the dotted and dot-dashed lines. The fact that the fixed-order behaviour at large q T is approximated better when Q is smaller is not unexpected. By varying Q, we smoothly set the transverse-momentum scale below which the resummed logarithmic terms are mostly effective; when Q is smaller, the resummation effects are confined to a range of smaller values of q T .
To quantify the resummation-scale uncertainty on the cross section at small and intermediate values of q T , we proceed as in the case of the renormalization and factorization scales, and we vary Q by a factor of 2 up and down from a reference value. We choose the reference value Q = M H /2, because of the better quality of the behaviour of the corresponding cross section at large q T . From Fig. 8, we see that at NLL+LO accuracy a scale variation between 1/4M H and M H produces a variation of the cross section of about ±15% in the region around the peak. At NNLL+NLO accuracy (Fig. 9) the resummation-scale dependence is much reduced: when Q varies between M H /4 and M H the change in the cross section at the peak is about ±5%, i.e., smaller than the corresponding uncertainty from variations of the renormalization and factorization scales (see Fig. 2).
Throughout this section we used the MRST2004 set [70] of parton distribution functions at NLO and NNLO. The NLO and NNLO parton densities from Alekhin are currently being updated [71]. The CTEQ [72] and GRV [73] groups do not include sets of NNLO parton densities. The parton distribution sets of MRST, Alekhin and CTEQ include estimates of experimental uncertainties, which lead to effects below to about 5% on the total cross section for Higgs boson production at the LHC. We do not expect significantly different results in the case of the q T cross section at the LHC, and we refer to Ref. [17] for results and discussions about the effects of available parton densities on the total cross section.
The quantitative predictions presented up to now are obtained in a purely perturbative framework. It is known (see e.g. Ref. [29] and references therein) that the transverse-momentum distribution is affected by non-perturbative (NP) effects, which become important as q T becomes small. In impact parameter space, these effects are associated to the large-b region. In our perturbative study the integral over the impact parameter turns out to be dominated by the region where b ∼ < 0.1-0.2 GeV −1 , larger values of b being strongly suppressed by the resummation of the logarithmic terms in the gluon form factor. Thus we do not expect particularly-large NP effects in the case of Higgs boson production at the LHC. This expectation is in agreement with the findings in Refs. [40][41][42][43][44].
A customary way of modelling NP effects in the case of DY lepton-pair production is to introduce an NP transverse-momentum smearing of the distribution. This is implemented by multiplying the b-space perturbative form factor by an NP form factor. Several different parametrizations of the NP form factor are available in the literature [63,[74][75][76][77]; the corresponding NP parameters are obtained from global fits to DY data.
In the case of Higgs boson production, the estimate of NP effects is obviously more uncertain, since we cannot exploit available experimental data. In Ref. [78] we studied the impact of NP contributions on the q T spectrum of the Higgs boson, by applying the DY NP corrections of Refs. [74][75][76] to our resummed results at NLL accuracy. We also considered the effect of rescaling the DY NP coefficients by the factor C A /C F , to take into account the different colour charges of the initial-state partons (qq in the DY process, gg in Higgs boson production) in the hard-scattering subprocess. Alternatively, we used the NP coefficients extracted in Ref. [43] from a fit of data on Υ production, a production process that is more sensitive to the gluon content of the colliding hadrons. All these different quantitative implementations of NP corrections, although certainly not fully justified, can give an idea of the size of the NP effects on the Higgs boson spectrum.
The results of Ref. [78] show that the impact of the NP effects on the NLL resummed distribution is definitely below 10% for q T ∼ > 10 GeV, and it decreases very rapidly as q T increases. Moreover, when q T ∼ < 10 GeV, different parametrizations of the NP terms can lead to sizeably different relative effects, as a consequence of our present ignorance on the absolute value of the NP contributions.
In view of these results, in the present paper we limit ourselves to considering a simple parametrization of the NP contributions. We multiply the b-space resummed component W H ab (b, M,ŝ) on the right-hand side of Eq. (10) by a NP factor, S N P , which includes a gaussian smearing of the form The NP coefficient g N P is varied in the range suggested by the study of Ref. [43]: g N P = 1.67-5.64 GeV 2 . Note that this procedure, with these values of g N P , well approximates the quantitative spread of NP effects found in Ref. [78] at NLL accuracy. In Fig. 11 we plot the effect of the NP smearing on our best perturbative predictions, as given by the results at NNLL+NLO accuracy. The inner plot shows the relative deviation from the NNLL+NLO perturbative result, as defined by the ratio where dσ N P N N LL+N LO is the NNLL+NLO cross section, dσ N N LL+N LO , supplemented with the NP form factor. We see that the NP effects give deviations from the purely perturbative result that are below 10% for q T ∼ > 5 GeV. Comparing the inset plots in Figs. 5 and 11, we also notice that the inclusion of higher-order contributions (going from NLL+LO to NNLL+NLO) and of NP contributions have a qualitatively similar effect at intermediate and small values of transverse momenta: both contributions make the distribution harder. At the quantitative level, ∆ is much smaller than K −1 when q T ∼ > 10 GeV, while ∆ and K −1 are comparable when q T ∼ < 10 GeV. This points towards a non-trivial interplay between higher-order perturbative effects and NP effects at fixed value of the Higgs boson mass.
In summary, the comparison of the NLL+LO and NNLL+NLO results from small (around the peak region) to intermediate (say, roughly, q T ∼ < M H /3) values of transverse momenta shows a Figure 11: The NNLL+NLO perturbative results supplemented with the NP form factor in Eq. (89). The upper (lower) curve at small q T is obtained with g N P = 1.67 GeV 2 (g N P = 5.64 GeV 2 ). nice convergence of the resummed QCD predictions for the q T spectrum of the Higgs boson at the LHC. From this comparison and from the effects of variations of the renormalization, factorization and resummation scales, we conclude that the perturbative QCD uncertainty of the NNLL+NLO results is uniformly of about 10% over this range of transverse momenta. The perturbative and NP uncertainty increases at smaller values of q T (see Figs. 5 and 11); the perturbative uncertainty increases also at larger values of q T [20][21][22]. The perturbative uncertainty on the NNLO cross section [14], as estimated in the same manner (i.e. by comparing the NLO and NNLO results, and performing scale variations), is about 15% [17]. Our results on the q T spectrum are thus fully consistent with those on the total cross section, since the bulk of the events is concentrated at small and intermediate values of the Higgs boson q T .

Conclusions
In this paper we have considered the transverse-momentum spectrum of generic systems of highmass M produced in hadron-hadron collisions. Following our previous work on the subject [33,1], we have illustrated and discussed in detail a perturbative QCD formalism that allows us to resum the large logarithmic contributions in the small-q T region (q T ≪ M) and to consistently match the ensuing result to the fixed-order contributions in the large-q T region (q T ∼ M). The main features of our approach, that make it different from other implementations of b-space resummation presented in the literature, are summarized below.
• The resummation is performed at the level of the partonic cross section. The parton distributions are thus evaluated at the factorization scale µ F , which has to be chosen of the order of the hard scale M. The resummation formula is then organized in a form that is in close analogy with the case of event shapes variables in hard-scattering processes [54][55][56][57] and threshold resummation in hadronic collisions [58,59]: the various classes of logarithmic contributions are controlled by the QCD coupling α S (µ 2 R ) evaluated at the renormalization scale µ R . This procedure naturally allows us to perform a systematic study of renormalizationand factorization-scale dependence, as is customarily done in fixed-order calculations. This should be compared with the other implementations of b-space resummation, where the scale at which the parton distributions are evaluated is of the order of 1/b, which also necessarily requires an extrapolation of the parton distributions in the NP region.
• The large logarithmic contributions are exponentiated in the form factor exp{G N }, where the function G N (see Eq. (14)) is universal: it does not depend on the produced high-mass system, and it only depends on the flavour of the partons involved in the hard-scattering subprocess. More precisely (see Appendix A), various process-independent form factors control the various partonic channels. The process dependence, as well as the factorizationscale and factorization-scheme dependence is fully included in the hard-scattering coefficient H N (see Eq. (12)).
• We impose a constraint of perturbative unitarity through the replacement in Eq. (16): the b-space form factor exp{G N ( L)} is equal to unity at b = 0. This constraint has a twofold purpose. On one hand, it avoids the introduction of unjustified higher-order contributions in the small-b region, which are present [79] in standard implementations of b-space resummation. On the other hand, it allows us to recover the total cross section at the nominal fixed-order accuracy upon integration over q T . Note that, as a consequence, perturbative uncertainties at intermediate values of q T are reduced.
The resummation formalism has been applied to the production of the SM Higgs boson in pp collisions. We combined the most advanced perturbative information that is available at present for this process: NNLL resummation at small q T and fixed-order perturbation theory at NLO at large q T . We developed a numerical code, named HqT [69], that performs the calculation at NLL+LO and NNLL+NLO accuracy. In Sect. 3.1 we have presented a selection of results that can be obtained by our program at LHC energies. Owing to the unitarity constraint, the integral of our spectra at NLL+LO (NNLL+NLO) correctly reproduces the total NLO (NNLO) cross sections. The results show a high stability with respect to scale variations and an increasing stability when going from NLL+LO to NNLL+NLO accuracy. As summarized at the end of Sect. 3.1, this suggests that the uncertainty from missing higher-order perturbative contributions is under good control.
A Appendix: Exponentiation in the multiflavour case In Sects. 2.2 and 2.3, we have discussed the exponentiation structure of the resummed component of the q T distribution. To simplify the notation and the presentation, we have limited ourselves to illustrating the case in which the partonic scattering involves a single flavour of partons. This appendix is devoted to generalize the exponentiation to the case with partons of different flavours.
To obtain Eq. (43), the multiflavour analogue of Eq. (12), we start from the representation in Eq. (42) of the resummed partonic cross section W F ab, N , and then we proceed as in Sect. 2.3. The main difference with respect to the steps in Eqs. (44)- (46) is that the solution of the QCD evolution equations (40) has the customary form † where the symbol P on the right-hand side denotes the path ordering expansion of the exponential matrix. Because of its matrix structure, the exponential in Eq. (91) has only a formal meaning. To recast Eq. (91) in a true exponential form, we can perform a systematic logarithmic expansion of the solution of the Altarelli-Parisi equations, by using a well-known method that dates back, at least, to Ref. [80].
The evolution operator in Eq. (91) can be written in the following form [80] (see also Ref. [81] for technical details): where U (LO) N is determined by the lowest-order anomalous dimensions γ and the operator V N fulfils the following differential equation: The evolution equation (93) can be solved by diagonalizing the anomalous dimensions matrix γ N , which has three different eigenvalues γ (1) i, N : one eigenvalue in the flavour non-singlet sector (i = NS), and two eigenvalues in the flavour singlet sector (i = ±). The solution of Eq. (93) is where E (i) N denotes the projector onto the flavour eigenspace corresponding to the eigenvalue γ (1) i, N . By inspection of Eq. (94), we see that it can be solved by performing a perturbative expansion, and the perturbative coefficients V (97) † In this appendix we use the boldface notation X to denote the flavour space matrix whose matrix elements are X ab = (X) ab .
We now come back to the right-hand side of Eq. (42). The evolution operator U N (b 2 0 /b 2 , µ 2 F ) is rewritten as U N (b 2 0 /b 2 , µ 2 F ) = U N (b 2 0 /b 2 , Q 2 ) U N (Q 2 , µ 2 F ). Then U N (b 2 0 /b 2 , Q 2 ) is replaced by the expression in Eq. (92). Equation (42) where we have defined the perturbative function and inside the curly brackets we have collected all the factors, S c , C N and U (LO) N , that depend on the impact parameter b. These factors contain the logarithmically-enhanced contributions that have to be resummed and organized in exponential form. The factor S c can be rewritten as where (see Eq. (31)) G c (α S (µ 2 R ), L; M 2 /µ 2 R , M 2 /Q 2 ) = − Q 2 b 2 0 /b 2 dq 2 q 2 A c (α S (q 2 )) ln M 2 q 2 + B c (α S (q 2 )) .
Note a key point: Eq. (104) does not regard the flavour matrix C N , but rather its matrix element C ca, N . Therefore, its right-hand side involves a true c-number exponential instead of a formal matrix exponential.
Inserting Eqs. (100), (102) and (104) From Eqs. The evolution operator U N (b 2 0 /b 2 , Q 2 ) does not contribute to the LL function g (1) (α S L) in Eq. (14). It starts to contribute to the resummation at the level of the NLL function g where U (NLO) N is the customary solution [80,81] of the evolution equations at NLO: The terms denoted by O(α n+2 S L n+1 ) on the right-hand side of Eq. (109) contribute at NNLL accuracy (they are of the same logarithmic accuracy as those in the function α S g (3) N (α S L) in Eq. (14)). Analogously, the terms denoted by O(α n+3 S L n+1 ) on the right-hand side of Eq. (109) contribute at NNNLL accuracy (they are of the same logarithmic accuracy as those in the function α 2 S g (4) N (α S L) in Eq. (14)). Therefore, to resum the NLL (NNLL) contributions to the form factor is sufficient to implement the solution of the evolution equations at the LO (NLO). Note, however, that, to be consistent with the resummed logarithmic expansion, the scale dependence of the running couplings α S (b 2 0 /b 2 ) and α S (Q 2 ) in Eq. (109) (Eq. (110)) has to be evaluated at the NLO (NNLO).

B Appendix: Bessel transformation of logarithmic contributions
where D(ǫ) = 2 b 0 2ǫ Γ(1 + ǫ) Note that the functions I n (x) exactly correspond to the following Bessel transformations: as can be checked by performing the limit q T → 0 of Eq. (112) or by verifying that the generating function in Eq. (127) has the following integral representation: The relation between I n (q T /Q) and the small-q T limit of I n (q T /Q) is not unexpected in view of the discussion in Sect. 2.2. The integral in Eq. (112) originates from Eq. (129) after the replacement L = ln(Q 2 b 2 /b 2 0 ) → L = ln(1 + Q 2 b 2 /b 2 0 ) at the integrand level: when q T → 0, such a replacement has no effects on the singular behaviour at any logarithmic accuracy.
Though I n (q T /Q) and I n (q T /Q) coincide when q T → 0, they behave quite differently at very large values of q T . When x → ∞, from Eqs. (116) and (127) Note, in particular, that I n (x) is integrable over x 2 when x → ∞, whereas I n (x) it is not.
The function I n (x) can easily be computed by performing the n-th derivative of the generating function (127) with respect to the parameter ǫ. To present the result, we first exclude the singular point x = 0 and consider only the region x > 0. Since the generating function depends on x only through the factor (x 2 ) −1−ǫ , x 2 I n (x) is simply a polynomial of degree n − 1 in the variable ln x 2 : where the coefficients d n are obtained from Eq. (128): The value of the first few coefficients is The choice x 0 = ∞ to define the plus-prescription in the case of I n is feasible since I n (x) (unlike I n (x)) is integrable over x 2 when x → ∞. This choice simplifies the definition of I n since the right-hand side of Eq. (142) (unlike Eq. (141)) does not contain any contact term proportional to δ(x 2 ). The contact term vanishes since the integrand factor ln n (1 + Q 2 b 2 /b 2 0 ) in Eq. (112) vanishes at b = 0. The vanishing of the contact term is thus ultimately related to the unitarity constraint in Eqs. (8) and (18).