Higgs-Boson Production at Small Transverse Momentum

Using methods from effective field theory, we have recently developed a novel, systematic framework for the calculation of the cross sections for electroweak gauge-boson production at small and very small transverse momentum q_T, in which large logarithms of the scale ratio m_V/q_T are resummed to all orders. This formalism is applied to the production of Higgs bosons in gluon fusion at the LHC. The production cross section receives logarithmically enhanced corrections from two sources: the running of the hard matching coefficient and the collinear factorization anomaly. The anomaly leads to the dynamical generation of a non-perturbative scale q_* ~ m_H e^{-const/\alpha_s(m_H)} ~ 8 GeV, which protects the process from receiving large long-distance hadronic contributions. We present detailed numerical predictions for the transverse-momentum spectrum of the Higgs boson, finding that it is quite insensitive to hadronic effects.


Introduction
The transverse-momentum spectra of electroweak bosons are among the most basic observables at hadron colliders. At large transverse momentum q T these spectra can be computed in fixedorder perturbation theory. On the other hand, if the transverse momentum is much smaller than the boson mass, higher-order corrections are enhanced by large Sudakov logarithms and must be resummed. This resummation was achieved a long time ago by Collins, Soper, and Sterman (CSS) [1] and has been implemented to high accuracy in several numerical codes [2][3][4][5][6]. In addition to the resummation at next-to-next-to-leading logarithmic (NNLL) accuracy, these programs include the matching to the next-to-next-to-leading (NNLO) fixed-order result [7]. For vector-boson production, we have revisited the resummation in the context of softcollinear effective theory (SCET) [8][9][10] and have derived an all-order factorization theorem for the cross section at small q T [11]. The factorization theorem for the differential cross section is affected by a collinear factorization anomaly, which generates an additional dependence on the hard momentum transfer in the low-energy theory. It was shown that in position space this extra dependence takes the form of a pure power of the gauge-boson mass. Relating our result to the traditional CSS formula then allowed us to determine the last missing ingredient needed for resummation to NNLL accuracy. The presence and all-order structure of these additional anomalous logarithms in the effective theory was confirmed by [12] and [13], but had been missed in earlier work on transverse-momentum resummation in SCET [14][15][16]. The work [12] derives the anomalous logarithms using a renormalization-group framework, in which the evolution is performed in rapidity instead of virtuality. Their final result agrees with the factorization formula derived in [11].
In the recent paper [17], we have used our formalism to analyze the differential cross section dσ/dq T at very small transverse momentum, where it exhibits quite remarkable properties. In this region, the spectrum is genuinely non-perturbative but dominated by short-distance physics and therefore calculable. For a boson of mass m V , long-distance effects are suppressed by a dynamically generated scale q * ∼ m V e −const/αs(m V ) , which is close to 2 GeV for the case of Z-boson production. While the underlying mechanism was identified a long time time ago [18], our framework has allowed us to systematically compute corrections also in this region. In [17] we have performed a detailed phenomenological study of Z-production at the Tevatron and the LHC and have investigated the numerical impact of long-distance effects.
The region of low transverse momentum is also important for the study of the Higgs boson and its properties, since the background is reduced when additional radiation is vetoed. In practice, this is done by imposing a jet veto. Several recent papers have considered resummation for the Higgs-boson cross section in the presence of such a cut [19][20][21][22][23], and the resummation is now known to NNLL accuracy [20,22]. In the present paper, we extend the formalism of [11] to the transverse-momentum spectrum of the Higgs boson and perform the resummation of the spectrum to NNLL accuracy. The Higgs case was also considered in [12], where expressions for the spectrum were presented at NLL order. Much of the extension is straightforward and some of the necessary ingredients were already given in [11], but there are a few interesting differences to the vector-boson case. First of all, while there is a single way to combine the spins of the incoming quarks to produce a vector boson, there are two ways to combine the gluon spins to produce a spin-zero state. As a consequence, the cross section is not just given by the product of two collinear functions as in the usual CSS formula, but a sum of two products of collinear functions describing the two production mechanisms, and the resummation formula must be modified accordingly [24]. We show that the collinear anomaly is the same for both structures, and that the dependence on the large scale m H therefore arises as an overall factor in position space. We then compute the collinear functions at one-loop order. The other feature, which distinguishes Higgs production from the vector-boson case, is that the infrared protection mechanism discussed above is much more efficient. The numerical value of q * ≈ 8 GeV is significantly higher than in the Z-boson case, and we show that long-distance hadronic effects have almost no impact on the Higgs-boson spectrum. We have implemented our resummed results for Drell-Yan, W , Z, and Higgs production in a public code CuTe [25] and give phenomenological predictions based on this program.

Factorization and resummation
We consider the cross section for the production of a Higgs boson with mass m H and transverse momentum q T = |q ⊥ | in gluon fusion at the LHC. The derivation of the factorization formula for the cross section proceeds exactly as in the case of the Higgs-production cross section defined with a jet veto, which we have recently considered in [20]. Our analysis there has been performed at fixed q ⊥ and rapidity y of the Higgs boson, and the integration over the boson phase-space was carried out at the end. We can thus immediately use the result for the factorized cross section obtained in [20], which reads where ξ 1,2 = √ τ e ±y and τ = (m 2 H + |q 2 ⊥ |)/s. The Born-level cross section is where √ s denotes the center-of-mass energy of the LHC and v is the Higgs vacuum expectation value. The Wilson coefficient C t multiplies the effective ggH operator obtained after integrating out the heavy top quark, while the hard matching coefficient C S arises when this operator is matched onto an effective two-gluon operator in SCET. Moreover, we have defined Here A c⊥ is the gauge-invariant effective gluon field of SCET, and S n , Sn denote soft Wilson lines. The soft function S describes the physics of soft gluons emitted from the colliding beam particles. The function B µν c is the standard transverse parton distribution function, first introduced in [26]. It describes the structure of the jet of collinear particles inside one of the colliding protons (the one moving along the light-like direction n µ ), which is probed at small transverse distance x ⊥ . The corresponding function B µν c for the second beam jet, consisting of anti-collinear particles (moving alongn), is given by the same formula with the replacements n → n and c →c. In [20], the sum over hadronic intermediate states was restricted by the jet veto, while in the present case the sums in (3) are completely inclusive.
In the context of SCET, proton matrix elements of collinear fields off the light cone are usually referred to as beam functions, a term introduced in [27] for fully unintegrated PDFs, which also depend on the other light-cone coordinate x − =n · x, in contrast to the matrix elements in (3). The two-loop renormalization of these functions was studied in [28], and the one-loop matching onto standard PDFs was calculated in [29]. The results of these two papers can, however, not be used in the present context, because in the limit x − → 0 light-cone singularities arise, which make the definition of the transverse PDFs in (3) subtle. In the context of SCET this was first discussed in [11], and we will now explain it in more detail.

Collinear anomaly
The beam-jet functions B µν c , B µν c and the soft function S suffer from light-cone divergences, which are not regularized by the conventional dimensional regularization procedure. These singularities cancel in the product of the three functions, but in order to make the individual objects well-defined an additional regularization is required [11]. At the end of the calculation the regulator can be removed, but a non-trivial effect remains: The additional regularization breaks a rescaling symmetry exhibited at the classical level by the effective Lagrangian of SCET (i.e., the property that the collinear matrix elements are invariant under a rescaling of the anti-collinear momenta and vice versa), which is not restored in the limit where the regulator is removed because it is spoiled by quantum effects. This "collinear factorization anomaly" manifests itself through a dependence of the product of soft and beam functions on the Higgs mass -the large scale in the problem. For the Drell-Yan process, the all-order form of this anomaly was derived in [11], and it was shown that the dependence on the hard scale (the electroweak gauge-boson mass in this case) takes the form of a pure power in x ⊥ space. We will now extend the derivation to the case of Higgs production and show that the form of the anomaly as well as the corresponding anomalous dimensions are spin independent.
The simplest way to introduce the additional regularization in such a way that gauge invariance and the factorization properties of the cross section are maintained is to regularize the phase-space integrals analytically by replacing [30] in the sum over intermediate states in (3). The scale ν is inserted to restore the canonical dimension of the integrals, in analogy to the scale µ of dimensional regularization. Instead of using a single light-cone component k + = k · n in the regulator term, one could also use the sum k + + k − = 2k 0 . This last form is similar to what is used in the "rapidity regularization scheme" proposed in [12]. For pertubative computations, using a single light-cone component is simplest, since SCET Feynman diagrams typically already contain light-cone denominators involving k + or k − . With the form (4) of the regularization, one finds that the soft function is given by scaleless integrals, and thus S(x ⊥ , µ) = 1 to all orders in perturbation theory [11]. We will no longer write it explicitly in the rest of the paper. If one were to use the energy k 0 instead of k + in (4), the soft function would be non-trivial, but can be absorbed into the beam-jet functions without loss of generality. The light-cone component k + of the anti-collinear particles moving along then µ direction is large, k + ∼ m H . Expanding in the regulator α, the dependence of the anti-collinear beam-jet function on the regulator scale ν thus takes the form ln(ν/m H ). On the other hand, the k + component of the collinear partons is small, In the collinear beam-jet function the dependence on the scale ν thus arises in the form ln(νx 2 T m H ). The requirement that the physical cross section (1) must be independent of the analytic regulator scale ν can then be expressed as In the factorization theorem (1) the Lorentz indices of the beam functions are contracted. The fact that ν independence also holds without contracting the indices follows by considering the factorization theorem for the production of a general color-neutral tensor field h µν . The corresponding factorization theorem has the same structure as (1), except that the hard matching coefficient |C S | 2 would now depend on the Lorentz indices of the tensor fields in the initial and final states. Since the logarithms in (5) have different arguments, the cancellation of the ν dependence among the different factors imposes a non-trivial constraint on the m H dependence of the product. As explained in detail in [11,31], the above equation implies that the dependence of the product of the two functions on m H must be power like. We can thus rewrite the product in the form The new beam-jet function B µν g (ξ, x ⊥ , µ) and the anomaly exponent F gg (x 2 T , µ) are independent of m H . Having determined the form of the anomaly, we now derive the scale dependences of the function B µν g (ξ 1 , x ⊥ , µ) and the exponent F gg (x 2 T , µ). Their anomalous dimensions can be inferred from the requirement that the cross section must be independent of the renormalization scale µ, which implies that the µ dependence of the product of beam functions must cancel against that of the hard function σ 0 C 2 t |C S | 2 in (1). This leads to the renormalization-group (RG) equations where Γ A cusp is the cusp anomalous dimension in the adjoint representation, and γ g is the anomalous dimension of the collinear gluon field as defined in [32]. The fact that each component of B µν g (ξ, x ⊥ , µ) renormalizes in the same way follows after considering the production of a general tensor field h µν and using that the anomalous dimension of the hard function is spin-independent [33]. Three-loop expressions for Γ A cusp and γ g can be found in [32]. Lorentz invariance implies that the renormalized beam-jet functions can be decomposed as where the coefficients only depend on the invariant The presence of the two different tensor structures (8) contributing to the Higgs cross section was first pointed out in [24]. As a result, the Higgs-production cross section does not simply factor into a product of two beam-jet functions, but instead it involves the above sum of products, which arises because the spin-0 state of two gluons is an entangled state. Rewriting the cross section (1) in terms of the new beam-jet functions, we obtain which is analogous to eq. (18) of our paper [11] on the Drell-Yan cross section for the production of electroweak gauge bosons. In the above formula, the three disparate scales m t ≫ m H ≫ x −1 T appear in factorized form, and it will be possible to resum logarithms of their ratios by controlling the µ dependence of the various functions in (10) using RG equations. After performing the Fourier integral the scales m H and q T get intertwined in a complicated way. At small q T this gives rise to interesting phenomena, such as a non-perturbative, but shortdistance dominated dependence of the cross section on q T /m H and α s (q T ) [11,17,34].

Refactorization
As long as the transverse displacement x T is much smaller than the scale of long-distance interactions in QCD (x T ≪ Λ −1 QCD ), formula (10) for the cross section can be simplified further, by matching the beam functions B (n) g onto ordinary parton distribution functions (PDFs), thereby computing their dependence on x 2 T in perturbation theory. The relevant matching relations read [1,26] which is valid up to hadronic corrections suppressed by powers of Λ 2 QCD x 2 T . The cross section can then be written in the final form whereC With a slight abuse of notation, we have traded the variables x 2 T and µ in the functions F gg and I without changing the names of these functions. This notation will be convenient for our discussion below, and it conforms with the notation used in [17]. Integrating the double differential cross section (12) over rapidity, we find where the parton luminosities and new kernel functions are defined as The factorized cross sections (12) and (15) receive power corrections in the two small quantities q 2 T /m 2 H and Λ 2 QCD x 2 T , which will not be indicated explicitly in our equations. A dependence on the hard scale m H enters formula (12) for the double-differential cross section in two places: via the hard matching coefficient C S and via an x T -dependent power of m H under the Fourier integral in (13). The latter effect is due to the collinear factorization anomaly [11]. As long as x 2 T ≪ Λ −2 QCD , the anomalous exponent F gg can be calculated in perturbation theory, and at least up to three-loop order it is related to the corresponding exponent F qq appearing in the Drell-Yan case by the Casimir-scaling relation [11] F gg (L ⊥ , a s ) Using the known expression for F qq , we then find where T F n f are the one-loop coefficients of the cusp anomalous dimension and β function, and contain the relevant two-loop information. Because of Casimir scaling, these coefficients take the same values as in the case of Drell-Yan production. At one-loop order, the kernel functions I (n) g←i (z, L ⊥ , a s ) are given by are the one-loop DGLAP splitting functions, and the remainder functions R g←i (z) and R ′ g←i (z) are given by The expression for I (1) g←i was calculated in [20], while the result for I (2) g←i is new. The oneloop functions I (i) g←g were also computed in Appendix C of [12]. 1 In [11], the SCET matching coefficients were related to the collinear functions in the traditional CSS formalism. Using these relations and the NNLO results presented in [7], it would be possible to extract the matching functions I (1) g←i (z, L ⊥ , a s ) at two-loop order.

Resummation of Sudakov logarithms
The resummation of large logarithms in the cross section (12) is accomplished by evolving the hard matching coefficients C t and C S to a scale µ at which the kernel functionsC gg←ij in (13) can be calculated using a controlled perturbative expansion. The solutions of the corresponding RG equations are discussed in detail in [20]. For our numerical work we use relations (A.2) and (A.4) from this paper. For the Drell-Yan case, the proper choice of the factorization scale µ has been discussed in [17]. Since the resulting expressions for Higgs production are completely analogous, we will not repeat details of the derivations here but rather summarize the main physical insights and quote the final expressions. The most naive choice would be to set µ ∼ x −1 T inside the Fourier integral in (13), in which case L ⊥ would be a small logarithm for any choice of x T . There are several disadvantages to such a treatment. First, since x T is integrated over all possible values, there would be no clear meaning to the scale µ in terms of a characteristic scale of the process. Second, setting the scale under the integral means that the integration unavoidably hits the Landau pole of the running coupling, giving rise to ambiguities in the numerical results. In the spirit of effective field theory, the scale µ should correspond to a physical scale in the underlying factorization theorem. We will choose it in such a way that on average the x T -dependent logarithm L ⊥ is small, and denote the corresponding value by µ ∼ x −1 T . Naively, one would expect that the transverse momentum q T and average transverse separation x T are conjugate variables satisfying q T ∼ x −1 T . While this is sometimes true, the general situation turns out to be more complicated. After integration over x ⊥ , the factorized dependence on m H and q T in (13) gets intertwined in a complicated way, and this gives rise to the peculiar effect that the two scales q T and µ ∼ x −1 T decouple for very small q T [17]. When this happens depends on the value of the coefficient As long as 0 < η < 1, one indeed finds that µ ∼ x −1 T ∼ q T , because contributions from large values x T ≫ q −1 T are suppressed due to the rapid oscillations of the phase factor of the Fourier integral, while contributions from small values x T ≪ q −1 T are phase-space suppressed. The situation changes, however, at very small transverse momentum, where with the prescription µ ∼ q T the value of η reaches 1. We denote by q * the value of µ where this happens, i.e.
where in the last step we have used the one-loop approximation for the running coupling. As long as q * is in the perturbative domain, one finds that at this scale x −1 T decouples from q T , and it remains a short-distance scale even in the extreme case where q T is taken to 0 [17]. Changing variables from x T to L ⊥ in the Fourier integral, one observes that the integrand exhibits a Gaussian peak with a width proportional to 1/ √ a s . The condition that at the peak , indicating that the factorization scale must be chosen in the vicinity of q * . We thus conclude that the proper scale choice is In our numerical work below, we will use µ = q T + q * as the default choice for the factorization scale. Solving the first equation in (24) numerically, we obtain q * ≈ 7.7 GeV for m H = 125 GeV, which is a short-distance scale well inside the perturbative domain. Due to the difference in color factors, this scale is significantly larger than in the case of the Drell-Yan production of electroweak gauge bosons, for which q * ≈ 1.75 GeV [17]. It follows from these arguments that the transverse-momentum distribution of Higgs bosons is protected from long-distance physics even for arbitrarily small q T -a fact that in the context of the Drell-Yan process has been pointed out first a very long time ago in [18]. The resummed perturbative series for the cross section generates the scale q * dynamically, and even though this is a short-distance scale, it is related to the boson mass m H in a genuinely non-perturbative way. The scale q * also sets the magnitude of hadronic long-distance corrections, which turn out to be power-suppressed in the ratio Λ QCD /q * . The dynamical origin of this suppression was studied in detail in [17]. We expect that these corrections are significantly smaller for Higgs production than for the Drell-Yan production of Z and W bosons. This expectation will be confirmed by our numerical studies presented below.
The above discussion shows that we must distinguish two regions of transverse momenta. For q T ≫ q * , the scale choice µ ∼ q T prevents that the logarithms L ⊥ give rise to large perturbative corrections. It is then consistent to count these logarithms as L ⊥ ∼ 1 and construct the perturbative series as a series in powers of a s . A different situation is encountered for q T ≪ q * . Even though the scale choice µ ∼ q * ensures that L ⊥ = O(1) on average, the Gaussian weight factor allows for significant contributions to the Fourier integral over a range of larger L ⊥ values with a width proportional to 1/ √ a s . It is then necessary to reorganize the perturbative expansion by adopting the modified power counting L ⊥ ∼ 1/ √ a s . This implies that single-logarithmic terms (a s L ⊥ ) n ∼ a n/2 s are always suppressed, whereas doublelogarithmic terms (a s L 2 ⊥ ) n ∼ 1 are unsuppressed and must be resummed to all orders. To keep track of this fact, we introduce an auxiliary expansion parameter ǫ (which at the end is set to 1) and assign the power counting a s ∼ ǫ and L ⊥ ∼ ǫ −1/2 . The terms contributing up to O(ǫ) to the cross section have been derived in [17] using recursive solutions of the relevant RG equations. Adapting the resulting expression to the present case, we find that the hard-scattering kernels defined in (13) can be written in the form where Note that we treat ln(m 2 H /µ 2 ) as a large logarithm and count η defined in (23) as an O(1) variable. The auxiliary parameter ǫ counts the order in a s resulting (for q T ≪ q * ) after the x T integral in (26) has been performed. The two terms given in the first line are unsuppressed and must be kept in the exponent of the integrand in (26), whereas the remaining terms can be expanded in powers of ǫ 1/2 . It is important that the expansion is truncated at an integer power of ǫ. The resulting integrals over the Bessel function can readily be evaluated numerically.
We finally give the expressions for the collinear kernel functions corresponding to our modified power counting. We find whileĪ (2) g←i coincides with I (2) g←i in (20) up to higher-order terms in ǫ. The corresponding contribution to (26) is of O(ǫ 2 ) and can be neglected to the order we are working. The quantities involve the convolutions of two DGLAP splitting functions. Following [17], we find (30) Equations (26)-(28) are our main results. With the help of these expressions, large logarithms can be resummed at NNLL order for arbitrarily small transverse momenta. For larger q T values, the additional terms contained in (27) and (28) reduce to higher-order terms proportional to a 2 s and a 3 s , which can be kept without doing any harm. Our formula thus provides a smooth interpolation between the regions of small and very small q T . In fact, it has been shown in [17] that the additional terms needed at very low q T serve an important purpose also at q T > q * , as they resum an asymptotically divergent, but Borel-summable subset of higher-order corrections in a s . In Figure 1, we present predictions for the resummed transverse-momentum distributions of Higgs bosons with mass m H = 125 GeV produced in gluon fusion at the LHC. The blue bands correspond to the results obtained at NLL order, in which case we keep terms up to O(ǫ 0 ) in (27) and (28) and use leading-order (LO) expressions for the Wilson coefficients C t and C S in RG-improved perturbation theory. The green bands show results with NNLL accuracy, which are obtained by retaining all terms shown in (27) and (28) and using next-to-leading-order (NLO) expressions for C t and C S . We use v = 246.675 GeV for the Higgs vacuum expectation value and work with MSTW2008NNLO PDFs with α s (M Z ) = 0.1171 [36]. In our numerical results, we include finite top-quark mass effects by working with the exact leading-order cross section, which can be found e.g. in [37]. Numerically, using m t = 172.6 GeV, this leads to a 6.6% increase in the cross section compared to the m t → ∞ limit. We do not include the b-quark contribution and electroweak effects, which are of similar size but opposite sign and tend to cancel each other. The combined effect of these would be a reduction of the cross section by about 1.5%. The left plot in the figure shows the spectrum dσ/dq T , while the plot on the right shows the distribution dσ/dq 2 T , which tends to a constant at the origin. Note that at NNLL order the value of the intercept is predicted with very good accuracy. In [17], we have presented an explicit formula for the intercept for the case of the Drell-Yan process.
To estimate the theoretical uncertainty of our predictions, we vary the factorization scale by a factor 2 about the default value µ = q T + q * . Our results also depend on the two hard matching scales µ t and µ h . At the scale µ t the top quark is integrated out, giving rise to an effective ggH vertex. The associated matching corrections are small, and varying µ t in the NNLL order cross section by a factor 2 about the default value µ t = m t changes the result by less than 1%. The scale µ h is associated with the hard momentum transfer in the ggH coupling, as described by the time-like gluon form factor. In [37,38], we have shown that it is advantageous to evaluate the relevant matching coefficient C S (−m 2 H , µ h ) with a time-like scale choice µ 2 h = −m 2 H . This eliminates the large perturbative corrections arising when the gluon form factor is continued from space-like to time-like kinematics. When this is done, also the corrections to C S (−m 2 H , µ h ) are of moderate size, and the effect of a variation of µ h by a factor 2 on the cross section is again below 1%. Since the variation of the factorization scale µ leads to the largest scale uncertainties by far, we use it to generate the error bands in the plots, keeping the hard matching scales µ t and µ h fixed at their default values. The NNLL corrections have a noticable effect and strongly enhance the cross section in the peak region. From the right plot, we observe that the scale variation at NLL order is very small in the vicinity of q T = 5 GeV, because the predictions with high and low µ values cross each other. Near such a band crossing, the scale variation underestimates the theoretical uncertainty, and it is therefore not too surprising that the NLL and NNLL bands do not overlap in the lowq T region. Since the prediction at NNLL order does not exhibit a band crossing, we believe that its scale variation provides a more reliable estimate of the theoretical uncertainty. Since we observe that the one-loop correction arising at NNLL order is larger than the NLL scale dependence suggests, we will be conservative when performing the matching to the fixed-order result in Section 3 and adopt a matching scheme which yields larger scale uncertainties for the combined result than for the resummed result itself (see Figure 4).
We proceed to study the impact on long-distance hadronic effects on the transversemomentum distribution, which for the case of Drell-Yan production of electroweak bosons are known to have a non-negligible impact [2,4,39]. Following [17], we model these effects by noting that the beam-jet functions B (10), which are nothing but transverseposition dependent PDFs, must vanish rapidly when the two gluon fields are separated by a transverse distance x T larger than the proton size. This motivates an ansatz of the form where the perturbative functions B (n) pert g carry all the scale dependence and are given by (11), whereas the hadronic form factor f hadr (r) with f hadr (0) = 1 describes the fall-off at large transverse distances and is parameterized in terms of a hadronic scale Λ NP . For simplicity, we will assume that this form factor is independent of ξ. The above ansatz inserts a factor [f hadr (x T Λ NP )] 2 under the integral over x T in (26), which suppresses the region of very large x T values. We will employ the Gaussian model for the form factor. For the case of Drell-Yan production, it was shown in [17] that the functional form of the model function only has a minor impact on the results, which are mainly sensitive to the value of the parameter Λ NP , and we have confirmed that the same is true in the present case. Choosing Λ NP ≈ 600 MeV shifts the position of the peak of the q T distribution for Z-boson production at the LHC from 3.2 GeV to 3.5 GeV and yields to a significantly better agreement with the data. A similar effect is seen for Tevatron data.
In Figure 2, we compare the situation in Drell-Yan production of Z bosons, for which the characteristic scale q * ≈ 1.75 GeV is rather low, with that in Higgs production at the LHC, for which q * ≈ 7.7 GeV is safely in the perturbative domain. As expected, we find that the impact of hadronic effects is significantly reduced in the latter case. With Λ NP ≈ 600 MeV, for instance, the peak position shifts by merely 100 MeV (from 9.1 GeV to 9.2 GeV), which is hardly visible on the scale of the plot. Even for Λ NP = 1 GeV, the shift amounts to only 300 MeV. We will see in the next section that perturbative uncertainties are significantly larger than this effect. As long as Λ NP is in the GeV range, it therefore seems safe to ignore the potential impact of long-distance effects for all practical purposes.

Predictions for the LHC
Having discussed the factorization of the cross section and its behavior at very small q T , we now present our final results for the transverse-momentum spectrum of Higgs bosons produced in gluon fusion at the LHC. In order to obtain reliable predictions also at intermediate q T values, we match the resummed differential cross section (15) to the O(α s ) fixed-order result. To this end, we add the fixed-order cross section σ NLO to the resummed result and subtract the fixed-order expansion of (15) so as to avoid double counting: The expansion of the resummed result to O(α s ) can be derived using g←j (z 2 ) .  (34), because they are to a large extent universal. This is obvious for the prefactor C 2 t , but it should also be the case for the large Sudakov logarithms and other corrections encoded in |C S | 2 , which also appear in the total cross section [37,38]. Factoring out the hard function H(µ) before performing the matching leads to a Sudakov suppression of the matching correction at low q T and an enhancement at larger transverse momentum.
In Figure 3, we compare our results for the resummed and matched differential cross section obtained when the hard function is not factored out ("naive matching") with the those found when the right-hand side of (34) is multiplied by H(µ) ("alternative matching"). In the second scheme, we indeed observe the Sudakov suppression of the matching correction at low q T , but also a drastic reduction of the scale dependence as displayed by the width of the band. This is due to a cancellation of the µ dependence between the resummed result and the matching correction. While some reduction may be expected, such a strong effect is presumably in part accidental. Indeed, instead of using H(µ) with variable µ we may equally well factor out the hard function at the fixed default scale µ = q T + q * . Doing so has the same qualitative effect on the matching correction, but as shown in the third plot the strong cancellation of scale dependence is not observed. To be conservative, we adopt this last choice as our "default matching" prescription.
Since the matching correction for Higgs-boson production is several times larger than that for the Drell-Yan case (see e.g. [17]), it would be preferable to extend the matching to the fixed-order cross section to two-loop order. This requires some effort, but it is possible since the corresponding fixed-order result is known [40][41][42] and has been implemented in several public codes, e.g. MCFM [43] and HNNLO [44]. We note that the quark beam function I q→q (z, x 2 T , µ) has recently been computed to two-loop accuracy [45]. Once this result is extended to the gluon channel, all two-loop ingredients for the resummed expression (15) will be known, and the matching should then be extended to O(α 2 s ). Our final results for the resummed and matched differential cross sections for Higgs pro- duction at the LHC, for √ s = 8 TeV and 13 TeV, are shown in Figure 4. The shape of the two spectra is very similar, the main effect of the higher center-of-mass energy being an increase in the cross section by about a factor of 2. The scale uncertainty is around ±10% in the peak region and increases for larger q T , as indicated in the panels below the plots. Our results are fully compatible with the NNLL order predictions of [4] obtained in the traditional CSS resummation framework [1]. The uncertainties found in that paper are slightly smaller in the peak region, but about a factor of 2 smaller at large q T . The reason for the reduced scale uncertainty is that this work implements matching to O(α 2 s ) as well as the hard-collinear two-loop corrections, which were calculated in [7].
We have implemented the results of our calculations into a public numerical code CuTe [25], which produces resummed and matched NNLL+NLO results not only for Higgs production, but also for the Drell-Yan process and the production of Z and W bosons at hadron colliders. The CuTe program computes cross sections with scale uncertainties and also includes different models for non-perturbative effects. It uses the LHAPDF interface, so it is straightforward to switch between different PDF sets. As an example, we show in Figure 5 the result obtained with NNPDF 2.3 PDF sets [46] and compare it with our default choice. The width of the bands now reflects the PDF uncertainties at the level of one standard deviation. We find good agreement of the distributions obtained with NNPDF 2.3 and MSTW2008 NNLO, the former giving rise to slightly higher cross sections. The bands correspond to one standard deviation variations of the PDF sets as compiled in MSTW2008 NNLO (green) and NNPDF 2.3 (orange). They do not include the uncertainty associated with the value of α s . For NNPDF, we choose the PDF set with the MSTW default value α s (M Z ) = 0.1170.

Conclusions
We have extended our previously developed formalism for the calculation of Drell-Yan cross sections at low and very low transverse momentum q T [17] to case of Higgs-boson production in gluon fusion at a hadron collider. Large Sudakov double logarithms in m H /q T are resummed in a systematic way using RG and anomaly equations. The leading logarithms are contained in the hard matching coefficient C S , which arises when the effective ggH operator is matched onto a two-gluon current operator in SCET. They are resummed by evolving this coefficient from a hard matching scale µ 2 h ∼ −m 2 H down to a characteristic scale µ ∼ q T . At NLL order, additional large logarithms are generated by the collinear factorization anomaly [11]. In transverse-position space, they give rise to an extra factor of the form (x 2 T m 2 H ) −Fgg(x 2 T ,µ) . Long-distance hadronic effects are controlled by a dynamically generated scale q * , which sets an upper limit on the average transverse separation in the process, x T 1/q * . As long as the value of q * lies in the perturbative domain, the transverse-momentum distribution is shortdistance dominated even for very small values of q T . Due to the difference in color factors (C A vs. C F ), this scale is significantly larger in the case of Higgs production (q * ≈ 7.7 GeV) compared with that of Z-boson production (q * ≈ 1.75 GeV). As a result, we find that longdistance hadronic corrections to the shape of the transverse-momentum distribution of Higgs bosons produced at the LHC have a negligible effect, see Figure 2.
From a technical point of view, a novel feature of the case of Higgs production is that the tensor structure of the beam-jet functions B µν g in (8) implies that the factorized cross sections (10) and (12) cannot be expressed as a simple product of two beam functions, but rather as a sum over products, reflecting the entanglement of the spin-0 state of two gluons [24]. However, we have shown that one can find a basis in which this effect arises first at NNLO in α s .
Our numerical analysis of the transverse-momentum distribution of Higgs bosons produced in gluon fusion at the LHC (with √ s = 8 and 13 TeV) presented in Section 3 gives rise to robust predictions, which hopefully will soon be confronted with first data from the LHC. We have explored different schemes for performing the matching to fixed-order perturbation theory and found that the various results agree within the estimated theoretical uncertainties. Our final predictions are presented in Figure 4, with PDF uncertainties estimated in Figure 5.
They suggest that at NNLL+NLO our predictions for the differential cross section dσ/dq T have uncertainties of order ±10% in the peak region, which seems reasonable given what is known about the behavior of the perturbative series for the total cross section. In order to reduce these uncertainties further, it would be necessary and desirable to extend our analysis to N 3 LL+NNLO, or at least to perform the matching to fixed-order theory at two-loop order.