Higgs boson production at the LHC: transverse-momentum resummation and rapidity dependence

We consider Higgs boson production by gluon fusion in hadron collisions. We study the doubly-differential transverse-momentum (q_T) and rapidity (y) distribution of the Higgs boson in perturbative QCD. In the region of small q_T (q_T<<M_H, M_H being the mass of the Higgs boson), we include the effect of logarithmically-enhanced contributions due to multiparton radiation to all perturbative orders. We use the impact parameter and double Mellin moments to implement and factorize the multiparton kinematics constraint of transverse- and longitudinal-momentum conservation. The logarithmic terms are then systematically resummed in exponential form. At small q_T, we perform the all-order resummation of large logarithms up to next-to-next-to-leading logarithmic accuracy, while at large q_T (q_T ~ M_H), we apply a matching procedure that recovers the fixed-order perturbation theory up to next-to-leading order. We present quantitative results for the differential cross section in q_T and y at the LHC, and we comment on the comparison with the q_T cross section integrated over y.


Introduction
The search for the Higgs boson [1] and the study of its properties (mass, couplings, decay widths) at hadron colliders require a detailed understanding of its production mechanisms. This demands reliable computations of related quantities, such as production cross sections and the associated distributions in rapidity and transverse momentum. In this paper we consider the production of the Standard Model (SM) Higgs boson by the gluon fusion mechanism.
The gluon fusion process gg → H, through a heavy-quark (mainly, top-quark) loop, is the main production mechanism of the SM Higgs boson H at hadron colliders. When combined with the decay channels H → γγ and H → ZZ, this production mechanism is one of the most important for Higgs boson searches and studies over the entire range, 100 GeV ∼ < M H ∼ < 1 TeV, of Higgs boson mass M H to be investigated at the LHC [2]. In the mass range 140 GeV ∼ < M H ∼ < 180 GeV, the gluon fusion process, followed by the decay H → W W → ℓ + ℓ − νν, can be exploited as main discovery channel at the LHC and also at the Tevatron [3], provided the background from tt production is suppressed by applying a veto cut on the transverse momenta of the jets accompanying the final-state leptons.
The dynamics of the gluon fusion mechanism in controlled by strong interactions. Detailed studies of the effect of QCD radiative corrections are thus necessary to obtain accurate theoretical predictions.
In QCD perturbation theory, the leading order (LO) contribution to the total cross section for Higgs boson production by gluon fusion is proportional to α 2 S , α S being the QCD coupling. The QCD radiative corrections to the total cross section are known at the next-to-leading order (NLO) [4]- [7] and at the next-to-next-to-leading order (NNLO) [8]- [12]. The Higgs boson rapidity distribution is also known at the NLO [13] and at the NNLO [14,15]. The effects of a jet veto have been studied up to the NNLO [11,14,15]. We recall that all the results at NNLO have been obtained by using the large-M t approximation, M t being the mass of the top quark. This approximation is justified by the fact that the bulk of the QCD radiative corrections to the total cross section is due to virtual and soft-gluon contributions [16,9,10,11,17]. The soft-gluon dominance also implies that higher-order perturbative contributions can reliably be estimated by applying resummation methods [9] of threshold logarithms, a type of logarithmically-enhanced terms due to multiple soft-gluon emission. In Ref. [17], the NNLO calculation of the total cross section is supplemented with threshold resummation at the next-to-next-to-leading logarithmic (NNLL) level; the residual perturbative uncertainty at the LHC is estimated to be at the level of better than ±10%. The NNLL+NNLO results [17] are nicely confirmed by the more recent computation [18]- [20] of the soft-gluon terms at N 3 LO; the quantitative effect [18] of the additional (i.e., beyond the NNLL order) single-logarithmic term at N 3 LO is consistent with the estimated uncertainty at NNLL+NNLO. The effect of threshold logarithms on the rapidity distribution of the Higgs boson has been considered in Ref. [21].
The gluon fusion mechanism at O(α 2 S ) produces a Higgs boson with a vanishing transverse momentum q T . A large (or, however, non-vanishing) value of q T can be obtained only starting from O(α 3 S ), when the Higgs boson is accompanied by at least one recoiling parton in the final state. This mismatch by a power of α S is a preliminary indication of the fact that the small-q T and large-q T regions are controlled by different dynamics regimes.
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. The LO, i.e. O(α 3 S ), calculation is reported in Ref. [22]. The results of Ref. [22] and the higher-order studies of Refs. [23,24] show that the large-M t approximation is sufficiently accurate also in the case of the q T distribution when q T ∼ < M H , provided q T ∼ < M t . Using the large-M t approximation, the NLO QCD computation of the q T distribution of the SM Higgs boson is presented in Refs. [25,26,27,14,15]. QCD corrections beyond the NLO are evaluated in Ref. [28], by implementing threshold resummation at the next-to-leading logarithmic (NLL) level. The results of the numerical programs of Refs. [14,15] can also be safely (i.e. without encountering infrared divergences) extended from large values of q T to q T = 0: in the small-q T region these programs evaluate the q T distribution up to NNLO.
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 definitely 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 ). The logarithmic terms are produced by multiple emission of soft and collinear partons (i.e. partons with low transverse momentum). 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 [29]- [36]. In the case of the SM Higgs boson, resummation has been explicitly worked out at leading logarithmic (LL), NLL [35,37] and NNLL [38] 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 at the LHC have been performed in Refs. [39]- [45], by combining resummed and fixed-order perturbation theory at different levels of theoretical accuracy. Other recent studies of various kinematical distributions of the SM Higgs boson at the LHC are presented in Refs. [46]- [50].
In Refs. [41,44] we studied the Higgs boson q T distribution integrated over the rapidity. In the small-q T region, the logarithmic terms were systematically resummed in exponential form by working in impact-parameter and Mellin-moment space. A constraint of perturbative unitarity was 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 . This constraint thus decreases the uncertainty in the matching procedure of the resummed and fixed-order contributions. Our best theoretical predictions were obtained by matching NNLL resummation at small q T and NLO perturbation theory at large q T . NNLL resummation includes the complete NNLO result at small q T , and the unitarity constraint assures that the total cross section at NNLO is recovered upon integration over q T of the transverse-momentum spectrum. Considering SM Higgs boson production at the LHC, we concluded [44] that the residual perturbative QCD uncertainty of the NNLL+NLO result is uniformly of about ±10% from small to intermediate values of transverse momenta.
In this paper we extend our study to include the dependence on the rapidity of the Higgs boson. Using the impact parameter and double Mellin moments, we can perform the extension by maintaining all the main features of the resummation formalism of Refs. [36,44]. We are then able to present results up to NNLL+NLO accuracy for the doubly-differential cross section in q T and rapidity at the LHC.
The paper is organized as follows. In Sect. 2 we recall the main aspects of the resummation formalism, and we illustrate the steps that are necessary to include the dependence on the rapidity in the q T resummed formulae. In Sect. 3 we apply the formalism to the production of the SM Higgs boson at the LHC, and we perform quantitative studies on the q T and rapidity dependence of the doubly-differential cross section. Some concluding remarks are presented in Sect. 4. Additional technical details on the double Mellin moments of the resummation formulae are given in Appendix A.

Rapidity dependence in q T resummation
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 Higgs boson H, accompanied by an arbitrary and undetected final state X. The centre-of-mass energy of the colliding hadrons is denoted by √ s. The rapidity, y, of the Higgs boson is defined in the centre-of-mass frame of the colliding hadrons, and the forward direction (y > 0) is identified by the direction of the momentum p 1 .
According to the QCD factorization theorem, the doubly-differential cross section for this process is 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σ ab are the partonic cross sections, 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 rapidity,ŷ, and the centre-of-mass energy,ŝ, of the partonic cross section (subprocess) are related to the corresponding hadronic variables y and s:ŷ with the kinematical boundary |ŷ| < ln ŝ/M 2 ( |y| < ln s/M 2 ) andŝ > M 2 (s > M 2 ).
The partonic cross section dσ ab is computable in QCD perturbation theory. Its power series expansion in α S contains the logarithmically-enhanced terms, (α n S /q 2 T ) ln m (M 2 H /q 2 T ), that we want to resum. To this purpose, we use the general (process-independent) strategy and the formalism described in detail in Ref. [44]. The only difference with respect to Ref. [44] is that the resummation is now performed at fixed values of the rapidity y, rather than after integration over the rapidity phase space. In the following we briefly recall the main steps of the resummation formalism, and we point out explicitly the differences with respect to Ref. [44].
We first rewrite (see Sect. 2.1 in Ref. [44]) the partonic cross section as the sum of two terms, The logarithmically-enhanced contributions are embodied in the 'resummed' component dσ a 1 a 2 is free of such contributions, and it can be computed by truncation of the perturbative series at a given fixed order (LO, NLO and so forth). In practice, after having evaluated dσ a 1 a 2 and its resummed component at a given perturbative order, the finite component dσ a 1 a 2 is obtained by the matching procedure described in Sects. 2.1 and 2.4 of Ref. [44].
The resummation procedure of the logarithmic terms has to be carried out [30]- [34] in the impact-parameter space, to correctly take into account the kinematics constraint of transversemomentum conservation. The resummed component of the partonic cross section is then obtained by performing the inverse Fourier (Bessel) transformation with respect to the impact parameter b. We write † dσ (res.) where J 0 (x) is the 0th-order Bessel function, and the factor W embodies the all-order dependence on the large logarithms ln(M H b) 2 at large b, which correspond to the q T -space terms ln(M 2 In the case of the q T cross section integrated over the rapidity, the resummation of the large logarithms is better expressed [36,44] by defining the N-moments W N of W with respect to z = M 2 H /ŝ at fixed M H . In the present case, where the rapidity is fixed, it is convenient (see e.g. Refs. [51,52]) to consider 'double' (N 1 , N 2 )-moments with respect to the two variables z 1 = e +ŷ M H / √ŝ and z 2 = e −ŷ M H / √ŝ at fixed M H (note that 0 < z i < 1). We thus introduce W (N 1 ,N 2 ) as follows: More generally, any function h(y; z) of the variables y (|y| < − ln √ z) and z (0 < z < 1) can be considered as a function of the two variables z 1 = e +y √ z and z 2 = e −y √ z. Thus, throughout the paper, the (N 1 , N 2 )-moments h (N 1 ,N 2 ) of the function h(y; z) are defined as Note that the double Mellin moments can also be obtained (see e.g. Ref. [53]) by introducing a Fourier transformation with respect to y (with conjugate variable ν = i(N 2 − N 1 )) and then performing a Mellin transformation with respect to z (with conjugate variable N = (N 1 + N 2 )/2): In the following equations, the functional dependence on the scales µ R and µ F is understood.
The convolution structure of the QCD factorization formula (2) is readily diagonalized by considering (N 1 , N 2 )-moments: where are the customary N-moments of the parton distributions. The use of Mellin moments also simplifies the resummation structure of the logarithmic terms in dσ can indeed be organized in exponential form as follows: where . is the Euler number) and, to simplify the notation, the dependence on the flavour indeces has been understood.
The structure of Eq. (10) is in close analogy to the cases of soft-gluon resummed calculations for hadronic event shapes in hard-scattering processes [54] and for threshold contributions to hadronic cross sections [55,51,56]. The function H (N 1 ,N 2 ) (which is process dependent) does not depend on the impact parameter b and, therefore, its evaluation does not require resummation of large logarithmic terms. It can be expanded in powers of α S as where σ 0 (α S , M H ) is the lowest-order partonic cross section for Higgs boson production. The form factor exp{G} is process independent ‡ ; it 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 → ∞. The functional dependence on b is expressed through the large logarithmic terms α n S L m with 1 ≤ m ≤ 2n. More importantly, all the logarithmic contributions to G with n + 2 ≤ m ≤ 2n are vanishing. Thus, the exponent G can systematically be expanded in powers of α S , at fixed value of λ = α S L, as follows: The term Lg (1) collects the leading logarithmic (LL) contributions α n S L n+1 ; the function g (2) resums the next-to-leading logarithmic (NLL) contributions α n S L n ; g (3) controls the next-to-next-toleading logarithmic (NNLL) terms α n S L n−1 , and so forth.
Note that we use the logarithmic variable L (see Eq. (11)) to parametrize and organize the resummation of the large logarithms ln(M H b) 2 . We recall the main motivations [44] for this choice. In the resummation region M H b ≫ 1, we have L ∼ ln(M H b) 2 and the use of the variable L is fully legitimate to arbitrary logarithmic accuracy. When M H b ≪ 1, we have L → 0 (whereas § ‡ More precisely, it depends only on the flavour of the colliding partons (see Appendix A).
§ As shown in Appendix B of Ref. [44] (see Eqs. (131) and (132) therein), after inverse Fourier transformation to q T space, the b-dependent functions ln n (M H b) 2 and L n lead to quite different behaviours at large q T . When ln(M H b) 2 → ∞ !) and exp{G(α S , L)} → 1. Therefore, the use of L reduces the effect produced by the resummed contributions in the small-b region (i.e., at large and intermediate values of q T ), where the large-b resummation approach is not justified. In particular, setting b = 0 (which corresponds to integrate over the entire q T range) we have exp{G(α S , L)} = 1: this property can be interpreted [44] as a constraint of perturbative unitarity on the total cross section; the dynamics of the all-order recoil effects, which are resummed in the form factor exp{G(α S , L)}, produces a smearing of the fixed-order q T distribution of the Higgs boson without affecting its total production rate.
The resummation formulae (10), (12) and (13) can be worked out at any given (and arbitrary) logarithmic accuracy since the functions H and G can explicitly be expressed (see Ref. [44]) in terms of few perturbatively-computable coefficients denoted by The key role of these coefficients to fully determine the structure of transverse-momentum resummation was first formalized by Collins, Soper and Sterman [34,32,36]. The present status of the calculation of these coefficients for Higgs boson production is recalled in Sect. 3.
In the case of the q T cross section integrated over the rapidity, Eq. (10) is still valid, provided the double (N 1 , N 2 )-moments are replaced by the corresponding single N-moments W N , H N , G N (see Sect. 2.2 in Ref. [44]). The relation between double and single moments can easily be understood by inspection of Eqs. (6)- (8). We see that setting ν = 0 in Eq. (8) is exactly equivalent to integrate the cross section over the rapidity. Therefore, the functions W N , H N , G N in Ref. [44] are obtained by simply setting Moreover, from the results presented in Ref. [44], we can straightforwardly obtain the functions H (N 1 ,N 2 ) and G (N 1 ,N 2 ) from the functions H N and G N . Roughly speaking, we simply have More precisely, these equalities are valid in the simplified case where there is a single species of partons (e.g. only gluons). In the following we comment on the physical picture that leads to Eq. (14). The generalization to considering more species of partons does not require any further conceptual steps: it just involves algebraic complications related to the treatment of the flavour indeces. The multiflavour case is briefly illustrated in Appendix A.
In the small-q T (large-b) region that we are considering, the kinematics of the Higgs boson is fully determined by the radiation of soft and collinear partons from the colliding partons (hadrons) in the initial state. The radiation of soft partons cannot affect the rapidity of the Higgs bosons. On the contrary, the radiation of partons that are collinear to p 1 (p 2 ), i.e. in the forward (backward) region, decreases (increases) the rapidity of the Higgs boson as a consequence of longitudinalmomentum conservation (see Eq. (3)). Since the emissions of collinear partons from p 1 and p 2 are dynamically uncorrelated (factorized from each other), correlations arise only from kinematics. The use of the (N 1 , N 2 )-moments exactly factorizes (see Eqs. (2) and (9)) the kinematical constraint of longitudinal-momentum conservation. It follows that the (N 1 , N 2 )-dependence of W (N 1 ,N 2 ) is given by the product of two functions (say, ) that depends only on N 1 or N 2 , respectively. If all the partons have the same flavour, the two functions should be equal, and Eq. (14) directly follows from [W (N 1 ,N 2 The formalism illustrated in this section defines a systematic 'order-by-order' (in extended sense) expansion [44] of Eq. (4): it can be used to obtain predictions with uniform perturbative accuracy from the small-q T region to the large-q T region. The various orders of this expansion are denoted ¶ as LL, NLL+LO, NNLL+NLO, etc., where the first label (LL, NLL, NNLL, . . . ) refers to the logarithmic accuracy at small q T and the second label (LO, NLO, . . . ) refers to the customary perturbative order at large q T . To be precise, the NLL+LO term of Eq. (4) is obtained by including the functions g (1) , g (2) and the coefficient H (1) (see Eqs. (13) and (12)) in the resummed component, and by expanding the finite (i.e. large-q T ) component up to its LO term. At NNLL+NLO accuracy, the resummed component includes also the function g N and the coefficient H (2) (see Eqs. (13) and (12)), while the finite component is expanded up to NLO. It is worthwhile noticing that the NNLL+NLO (NLL+LO) result includes the full NNLO (NLO) perturbative contribution in the small-q T region.
We recall [44] that, due to our actual definition of the logarithmic parameter L in Eq. (10) and to our matching procedure with the perturbative expansion at large q T , the integral over q T of the q T cross section exactly reproduces the customary fixed-order calculation of the total cross section. This feature is not affected by keeping the rapidity fixed. Therefore, the NNLO (NLO) result for total cross section at fixed y is exactly recovered upon integration over q T of the NNLL+NLO (NLL+LO) q T spectrum at fixed y.
Within our formalism, resummation is directly implemented, at fixed M H , in the space of the conjugate variables N 1 , N 2 and b. To obtain the cross section in Eq. (2), as function of the kinematical variables s, y and q T , we have to perform inverse integral transformations. These integrals are carried out numerically. We recall [44] that the resummed form factor (i.e., each of the functions g (k) (α S L) in Eq. (13)) is singular at the value of b where α S (µ 2 R ) L = π/β 0 (β 0 is the first-order coefficient of the QCD β function). This singularity has its origin from the presence of the Landau pole in the running of the QCD coupling α S (q 2 ) at low scales. When performing the inverse Fourier (Bessel) transformation with respect to the impact parameter b (see Eq. (5)), we deal with this singularity by using a 'minimal prescription' [56,57]: the singularity is avoided by deforming the integration contour in the complex b space (see Ref. [57]). We note that the position of the singularity is completely independent of the values of N 1 and N 2 . Thus, the inversion of the Mellin moments is performed in the customary way (in Mellin space there are no singularities for sufficiently-large values of Re N 1 and Re N 2 ). In this respect, going from single N-moments (as in Ref. [44]) to double (N 1 , N 2 )-moments (as in the present case, where the rapidity is kept fixed) is completely straightforward, with no additional (practical or conceptual) complications.

Higgs boson production at the LHC
In this section we apply the resummation formalism of Sect. 2 to the production of the Standard Model Higgs boson at the LHC. We closely follow our previous study of the single differential (with respect to q T ) cross section, with the same choice of parameters as stated in Sect. 3 of ¶ In the literature on q T resummation, other authors sometime use the same labels (NLL, NLO and so forth) with a meaning that is different from ours.
We recall that the LO term at small q T (i.e. including the region where q T = 0) is proportional to α 2 S , whereas the LO term at large q T is proportional to α 3 S . This mismatch of one power of α S (and the ensuing mismatch of notation) persists at each higher order (NLO, NNLO, ...).
Ref. [44]. Therefore, the integration over y of the double differential (with respect to y and q T ) cross sections presented in this section returns the q T cross sections of Ref. [44]. As a cross-check of the actual implementation of the calculation, we have verified that after integration over the rapidity the numerical results in Ref. [44] are reobtained within a high accuracy.
As in Refs. [17,44], we use an 'improved version' [16] of the large-M t approximation. The cross section is first computed by using the large-M t approximation. Then, it is rescaled by a Born level factor, such as to include the exact lowest-order dependence on the masses, M t and M b , of the top and bottom * * quarks, which circulates in the heavy-quark loop that couples to the Higgs boson. We use the values M t = 175 GeV and M b = 4.75 GeV. As discussed in Ref. [17] and recalled in Sect. 1, this version 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.
For the sake of brevity, we present quantitative results only at NNLL+NLO accuracy, which is the highest accuracy that can be achieved by using the present knowledge of exact perturbative QCD contributions (resummation coefficients and fixed-order calculations [25]- [27]). We use the MRST2004 set [58] [35] and H (1) [37]. To reach the NNLL+NLO accuracy, the form factor function G (N 1 ,N 2 ) in Eq. (13) must include the contribution from g (3) (N 1 ,N 2 ) (which is controlled by the coefficients B (2) [38] and A (3) [59]), and the coefficient function H (N 1 ,N 2 ) in Eq. (12) has to be evaluated up to its second-order term H (2) (N 1 ,N 2 ) . In Ref. [44] we exploited the unitarity constraint G(α S , L)| b=0 = 0 to numerically derive an approximated form of the coefficient H (2) from the NNLO calculation [12] of the total cross section. The recent calculation of Ref. [15], which is based on the complete evaluation of H (2) (N 1 ,N 2 ) in analytic form, allows us to gauge the quality of the approximated form. We find that the use of the H (2) of Ref. [44] leads to differences of about 1% with respect to the exact computation of the rapidity cross section at NNLO.
All the numerical results in this section are obtained by fixing the renormalization and factorization scales at the value µ R = µ F = M H . The 'resummation scale' Q (the auxiliary scale introduced in Ref. [44] to gauge the effect of yet uncalculated logarithmic terms at higher orders) is also fixed at the value Q = M H . The mass of the Higgs boson is set at the value M H = 125 GeV.
We start our presentation of the predictions for Higgs boson production at the LHC by considering the q T dependence of the cross section at fixed values of the rapidity. In Fig. 1, we set y = 0 and we compare the customary (when q T > 0) NLO calculation (dashed line) with the resummed NNLL+NLO calculation (solid line).
As expected, the NLO result diverges to −∞ as q T → 0 and, at small values of q T , it has an unphysical peak that is produced by the numerical compensation of negative leading logarithmic and positive subleading logarithmic contributions. The presence of this peak is not accidental. At * * We note that the Born level cross section is not insensitive to the contribution of the bottom quark. Adding  large q T , the perturbative expansion at any fixed order has no pathological behaviour: it leads to a positive cross section, whose value decreases as q T increases. When q T → 0, instead, any fixedorder calculation diverges alternatively to ±∞ depending on the perturbative order. Therefore, to go smoothly from the large-q T behaviour to the small-q T limit, the NLO (or N 3 LO, and so forth) calculation of the cross section has to show at least one peak in the intermediate-q T region.
We recall once more that the label NLO in Fig. 1 refers to (and originates from) the perturbative expansion at large q T . To avoid possible misunderstandings (coming from such a label) when interpreting the dashed (NLO) curve in the small-q T region, we point out that, the only difference produced in Fig. 1 by the NNLO calculation at small q T (this calculation can be carried out, for example, by using the NNLO codes of Refs. [14,15]) is a spike around the point q T = 0. More precisely, as long as q T = 0, the dashed curve is exactly the result of the NNLO calculation of the q T cross section at small q T . The only difference introduced in the plot by this NNLO calculation would occur in the first bin (with arbitrarily small size) that includes the point q T = 0. The NNLO value of the q T cross section in this first bin is positive and fixed by the value of the NNLO total cross section † † . Of course, owing to the increasingly negative behaviour of the q T distribution when q T → 0, the NNLO value of the q T cross section in the first bin increases by decreasing the size of that bin.
The resummed NNLL+NLO result in Fig. 1 is physically well-behaved at small q T (it vanishes † † By definition, the integral over q T of d 2 σ/(dq T dy) at NNLO is equal to dσ/dy at NNLO. as q T → 0 and has a kinematical peak at q T ∼ 12 GeV), and it converges to the expected NLO result only when q T is definitely large (q T ≃ M H ).
To quantify more clearly the effect of the resummation on the the NLO result, the value at y = 0 of the q T dependent K-factor, is shown in the inset plot of Fig. 1. The dashed line shows the analogous K-factor as computed from the ratio of the rapidity integrated cross sections. The similarity between these two K-factors is a first indication of the mild rapidity dependence of the resummation effects. By inspection of the inset plot, we note that NNLL resummation is relevant not only at small q T , but also in the intermediate-q T region: as soon as q T ∼ < 80 GeV, the resummation effects are larger than 20%. Of course, the fact that K ∼ 1 at q T ∼ 24 GeV is purely accidental: it simply follows from the unphysical behaviour of the fixed-order perturbative expansion at small q T . Considering other values of the rapidity, from the central to the off-central rapidity region, we find the same features as observed at y = 0. Our results of the q T spectrum at y = 2 are presented in Fig. 2. The NNLL+NLO spectrum has a peak at q T ∼ 11 GeV. As happens in the case of the q T distribution integrated over y, the effect of NNLL resummation is definitely non-negligible starting from relatively-high values of q T . For example, at q T = 50 GeV the NNLL+NLO result is about 30% higher than the NLO result.
To analyze the rapidity dependence in more detail, we study the doubly-differential cross section at fixed values of q T . In Figs. 3 and 4, we show quantitative results at two typical values of the transverse momentum, q T = 15 GeV and q T = 40 GeV, in the small-q T and intermediate-q T region, respectively.  Figure 3 shows the rapidity distribution at NNLL+NLO (solid line) and NLO (dashes) accuracy when q T = 15 GeV. At this value of q T , the effect of NNLL resummation reduces the cross section. For example, when y = 0 the reduction effect is about 25%. As can be observed in the inset plot, the relative contribution from the resummed logarithmic terms is rather constant in the central rapidity region, and its dependence on y only appears in forward (and backward) region, where the cross section is quite small.
When q T = 40 GeV (see Fig. 4), instead, the effect of NNLL resummation increases the absolute value of the cross section. For example, when y = 0 the NLO cross section is increased by about 22%. Nonetheless, as for the relative effect of resummation and the rapidity dependence of the K-factor, we observe features that are very similar to those in Fig. 3. The resummation effects have a very mild dependence on y in the central and (moderately) off-central regions, and this explains the remarkable similarity between the solid and dashed lines in the inset plot of Figs. 1 and 2. Since the kinematical region where |y| ∼ < 2 accounts for most of the total cross section, when comparing the ratio K(q T , y) to the analogous ratio of the y-integrated cross sections, hardly any differences are expected, unless the large-rapidity region is explored.
The mild rapidity dependence of the q T shape of the resummed results can be studied with a finer resolution by defining the following ratio: This ratio gives the doubly-differential cross section normalized to the q T cross section integrated over the full rapidity range. For comparison, we consider also the q T -integrated version of the cross section ratio in Eq. (16), and we define the ratio of the rapidity cross section dσ/dy over the total cross section σ.
We have computed the ratio in Eq. (16) by using the resummed q T cross sections at NNLL+NLO accuracy. The results, as a function of q T , are presented in Fig. 5 (solid lines) at two different values, y = 0 and y = 2, of the rapidity. The results of the analogous (q T -independent) ratio R y (computed ‡ ‡ at NNLO with the numerical programs of Refs. [14,15]) at the corresponding values of rapidity are also reported (dotted lines) in Fig. 5. The dashed lines in Fig. 5 correspond to the computation of Eq. (16) by using the q T cross sections at NLO: we see that the dashed and solid lines are very similar (as expected from the similarity of the dashed and solid lines in the inset plot of Figs. 1 and 2). As discussed below, the results in Fig. 5 show that the cross section decreases and the q T spectrum softens when the rapidity increases.
We observe that the lines at y = 0 lie above the lines at y = 2; this is just a consequence of the fact that the cross sections (both at fixed q T and after integration over q T ) decrease when y increases.
At fixed y, R(q T ; y) is not constant: it depends (though very slightly) on q T . We note that the corresponding upper and lower lines in Fig. 5 have different slopes with respect to q T : fixing q T , the q T slope of R(q T ; y) decreases from positive to negative values as y increases from y = 0 to y = 2, thus showing that the q T spectrum becomes slightly softer at larger rapidity.
In general, as |y| increases, the hardness of the q T shape of d 2 σ/(dq T dy) decreases. Since the cross section decreases by increasing the rapidity, the hardness of dσ/dq T (the denominator in Eq. (16)) is intermediate between the values of the hardness of d 2 σ/(dq T dy) (the numerator in Eq. (16)) at y = 0 and at large |y|. As a consequence, the q T slope of R(q T ; y) is necessarily positive when y = 0. Note that the q T slope is already negative when y = 2 (Fig. 5): this is a consequence of the fact that the bulk of the cross section is in the rapidity region |y| ∼ < 2.
Our qualitative illustration of the results in Fig. 5 can be accompanied by some quantitative observations. We note that the rapidity dependence of the cross sections is sizeable: going from y = 0 to y = 2, the ratio R y decreases by about 43%; comparable variations affect the ratio R(q T ; y), which is not very different from R y and it is slowly dependent on q T . Indeed, at fixed y, the ratio R(q T ; y) at NNLL+NLO accuracy has a small and nearly constant slope from low values of q T around the peak (say, q T ∼ 10 GeV) to q T = 100 GeV; varying q T in this region, R(q T ; y) increases by about 11% when y = 0, and it decreases by about 16% when y = 2. In the same range of q T and y, the values of R(q T ; y) at NNLL+NLO (solid lines) and at NLO (dashed lines) are very similar: although this is expected at large q T , the differences never exceed the level of about 4% even at values of q T as low as q T ∼ 10 GeV.
In summary, the results in Fig. 5 show that, when |y| increases from the central to the (moderately) off-central region, the cross sections vary more in absolute value than in q T shape. These features deserve some words of discussion.
We first consider the total cross section σ and the rapidity cross section dσ/dy. We recall (see Sect. 1) that the value of these cross sections is sizeably affected by QCD radiative corrections. The bulk of the effect is due to the radiation of virtual and soft gluons, and they cannot affect the rapidity of the Higgs boson. As a consequence, the ratio R y has little sensitivity to perturbative QCD corrections. The decreases of R y as |y| increases is mainly driven by the decrease of the gluon density f g (x, M 2 H ) as x increases. Considering the large-q T region, similar arguments apply to the q T cross sections dσ/dq T and dσ/(dq T dy), and similar conclusions apply to the ratio R(q T , y). In the small-q T region, we have to consider the additional and large effect produced on the q T cross sections by the logarithmically-enhanced terms ln m (M 2 H /q 2 T ). These terms are due to the radiation of soft and collinear partons. As already discussed in Sect. 2, the rapidity of the Higgs boson can be varied only by collinear radiation, while soft radiation can only lead to on overall (independent of y) rescaling of the q T cross sections. At the LL level, only soft radiation contributes (the LL function g (1) in Eq. (13) does not depend on N 1 and N 2 ) and all the logarithmic terms cancel in the ratio R(q T , y). The y sensitivity of R(q T , y) starts at the NLL level. The corrections produced on the dominant soft-gluon effects by the collinear radiation are physically [29] well approximated by varying the scale µ of the gluon density from µ ∼ M H to µ ∼ q T . As a consequence, the variations of the hardness of the of the q T cross sections are mainly driven by d ln f g (x, q 2 T )/d ln q 2 T , the amount of scaling violation of the gluon density. Since the scaling violation decreases as x increases, the hardness of dσ/(dq T dy) decreases and the q T spectrum softens as |y| increases. Note that, by increasing x, the gluon density decreases faster than its scaling violation: this explains why dσ/(dq T dy) varies more in absolute value than in q T shape when |y| increases.
We conclude this section with some comments about the theoretical uncertainties on the doubly-differential cross section dσ/(dq T dy) at NNLL+NLO accuracy. In Ref. [44] the perturbative QCD uncertainties on dσ/dq T were investigated by comparing the results at NNLL+NLO and NLL+LO accuracies and by performing scale variations at NNLL+NLO level. We also considered the inclusion of non-perturbative contributions, and we found that they lead to small corrections provided q T is not very small. From these studies we concluded that the NNLL+NLO result has a QCD uncertainty of about ±10% in the region from small (around the peak of the q T distribution) to intermediate (say, roughly, q T ∼ < M H /3) values of transverse momenta. Similar studies can be carried out in the case of the doubly-differential cross section dσ/(dq T dy). These studies are not reported here since their results and the ensuing conclusions are very similar to those in Ref. [44]. The reason for this similarity is a feature that we have pointed out throughout this section: the q T resummation effects have a very mild dependence on the rapidity and, thus, they are almost unchanged when comparing dσ/(dq T dy) with dσ/dq T (equivalently, they largely cancel in the ratio in the cross section ratio of Eq. (16)).

Summary
We have considered the resummation of the logarithmically-enhanced QCD contributions that appear at small transverse momenta when computing the q T spectrum of a Higgs boson produced in hadron collisions. In our previous work on the subject [41,44], the rapidity of the Higgs boson was integrated over: resummation was implemented by using a formalism based on a transform to impact parameter and Mellin moment space. In this paper we have extended the resummation formalism to the case in which the rapidity is kept fixed, and we have considered the doublydifferential cross section with respect to the transverse momentum and the rapidity. We have shown that this extension can be carried out without substantial complications: it is sufficient to enlarge the conjugate space by introducing a suitably-defined double Mellin transformation.
The main aspects of our method [36,44], which are recalled here, are unchanged by the inclusion of the rapidity dependence. The resummation is performed at the level of the partonic cross section, and the parton densities are factorized as in the customary fixed-order calculations. The formalism is completely general and it can be applied to other processes: the large logarithmic contributions are universal and, thus, they are systematically exponentiated in a process-independent form (see Eqs. (10) and (13)); the process-dependent part is factorized in the hard-scattering coefficient H. A constraint of perturbative unitarity is imposed on the resummed terms (see Eq. (11)), so that the q T smearing produced by the resummation does not change the total production rate. This constraint reduces the effect of unjustified higher-order contributions at intermediate q T and facilitates the matching procedure with the complete fixed-order calculations at large q T . In particular, when the rapidity is kept fixed, the integration over q T of dσ/(dq T dy) at NNLL+NLO accuracy returns dσ/dy at NNLO.
We have presented numerical results for Higgs boson production at the LHC. Comparing fixedorder and resummed calculations, we find that the resummation effects are large at small q T (as expected) and still sizeable at intermediate q T . The inclusion of the rapidity dependence has little quantitative impact on this picture since, as we have shown, the q T resummation effects are mildly dependent on the rapidity. Going from the central to the (moderately) off-central rapidity region, the q T shape of the spectrum slightly softens. In the range from small to intermediate values of q T , the residual perturbative uncertainty of the NNLL+NLO predictions for dσ/(dq T dy) is comparable to that of advanced (NNLO or NNLL+NNLO) calculations of the q T inclusive cross sections dσ/dy and σ.

A Appendix
In this appendix we present the structure of the resummation formula (10) by explicitly including the dependence on the flavour indeces of the colliding partons.
In the context of our resummation formalism, a detailed derivation of exponentiation in the multiflavour case is illustrated in Appendix A of Ref. [44]. Considering a generic LO partonic subprocess c +c → F (F = H and c =c = g in the specific case of Higgs boson production by gluon fusion), and performing q T resummation after integration over the rapidity, the resummed component dσ (res.) a 1 a 2 /dq 2 T of the partonic cross section is controlled by the N-moments W F a 1 a 2 , N . The final exponentiated result for these N-moments is given by the master formulae (106)-(108) of Ref. [44]. We recall the master formula (106) in the following form: where the sum extends over the following set of flavour indices: and, for simplicity, the functional dependence on various scales (such as the renormalization and factorization scales) is understood. The functions G {I}, N and H {I}, F a 1 a 2 , N are given in the master formulae (107) and (108), respectively.
In the present paper, q T resummation is performed at fixed values of the rapidity, and the double (N 1 , N 2 )-moments W (N 1 ,N 2 ) F a 1 a 2 in Eq. (6) replace the N-moments W F a 1 a 2 , N of Ref. [44]. The generalization of Eq. (10) to the multiflavour case is straightforwardly obtained from Eq. (18) by the simple replacement N → (N 1 , N 2 ): The expressions in Eqs. (21) and (22) are completely analogous to the master formulae (107) and (108) in Ref. [44] (the functional dependence on the scales M, µ R , µ F and Q is explicitly denoted in those formulae). In particular, we note that the dependence of G (N 1 ,N 2 ) and H (N 1 ,N 2 ) on the Mellin variables N 1 and N 2 is completely factorized: each of terms on the right-hand side of Eqs. (21) and (22) depends only on one Mellin variable (either N 1 or N 2 ). This factorized structure is completely consistent with Eq. (14) and with the physical picture discussed below Eq. (14); the dependence on N 1 (N 2 ) follows the longitudinal-momentum flow and the flavour flow a 1 → b 1 → i 1 → c (a 2 → b 2 → i 2 →c) that are produced by collinear radiation from the initial-state parton with momentum p 1 (p 2 ). The various Mellin functions (G i, N , E (i) N , U N and so forth) in Eqs. (21) and (22) can be found in Ref. [44].