Precise QCD Description of the Higgs Boson Transverse Momentum Spectrum

The transverse momentum ($p_T$) distribution of Higgs bosons at hadron colliders enables a detailed probe of its production dynamics and is a key ingredient to precision studies of Higgs boson properties, but receives very large QCD corrections. We obtain a precision prediction for the $p_T$ spectrum by matching second-order (NNLO) QCD corrections at large $p_T$ with resummation of third-order logarithmic (N$^3$LL) corrections at small $p_T$. We achieve significantly improved results for $p_T<35\,{\rm GeV}$ with perturbative uncertainties $\lesssim \pm 6\%$, and thus a convergent perturbative series for all values of $p_T$.

Introduction.-With the increased statistics and excellent performance of the experiments at the Large Hadron Collider (LHC), precision analysis of the production and decay of the Higgs boson [1][2][3][4], the linchpin of the Standard Model, is now becoming reality. Concurrent with the experimental advances has been a concerted effort to improve the theoretical calculations of Higgs production through the gluon fusion mechanism , driven by the large corrections one finds in the first few orders of perturbation theory. A key observable is the transverse momentum (p T ) spectrum of the Higgs boson with respect to the beam directions, which quantifies how the Higgs boson recoils against partonic radiation, thereby probing the interaction of the Higgs boson with Standard Model particles and possible new states. Calculating the full spectrum requires a combination of fixed order perturbation theory in α s and a resummation to all orders in α s of the most singular large logarithms ln(p T /m H ) for p T m H , where m H is the Higgs mass. In this letter, we compute the transverse momentum spectrum under the infinite top quark mass assumption at the highest precision currently possible, by combining next-to-next-toleading order (NNLO) at large transverse momentum [20] with an all-order resummation (based on soft-collinear effective field theory, SCET [29][30][31][32]) of large logarithmic corrections at small transverse momentum, including sub-leading logarithms up to the third level (N 3 LL) using the recently computed 3-loop rapidity anomalous dimension [33,34]. Our numerical NNLO results extend to considerably smaller values of p T compared to earlier work [17,20]. Resummation at N 3 LL was previously achieved in a momentum space resummation using the RadISH Monte Carlo program [27] and matched to NNLO results at low resolution [17]. Our approach to obtain N 3 LL resummation uses a SCET factorization formula, which yields analytic expressions for all singular fixedorder terms. Taken together, these two advancements enable a detailed cross-validation of the results from fixedorder and resummation in all parton-level channels, and yield the first high-resolution description of the Higgs boson p T spectrum.
Method.-The dominant Higgs boson production process at the LHC is gluon fusion, mediated by a top quark loop. This process can be described by integrating out the top quark, resulting in an effective field theory (EFT) coupling the Higgs boson to the gluon field strength tensor [35][36][37], combined with QCD containing five massless quark flavours. The Wilson coefficient of this effective operator is known to three loop accuracy [38], and the validity of this description up to p T of the order of the top quark mass has been established [21,23] in detail. For higher p T the exact top quark mass dependence needs to be included, and a recent NLO QCD calculation [24] has demonstrated that this can be accounted for by a multiplicative rescaling of the EFT predictions. The results presented here focus on p T below the top quark threshold, and are obtained in the EFT framework.
The p T distribution of the Higgs boson, valid at both small and large p T , can be written as where dσ r /dp 2 T is the resummed distribution, dσ f /dp 2 T is the fixed order distribution, and dσ s /dp 2 T is the fixedorder singular distribution. The resummed distribution dσ r /dp 2 T can be written as an inverse Fourier transformation from impact parameter space b to momentum space p T , as given in Eq. (2). Here, σ 0 is the Born cross section for gg → H and we denote | b| = b and b 0 = 2e −γ E . dσ r dp 2 In Eq. (3), the kernel of the integral, W , has been factorized into products of a hard function H(m H , µ h ), a soft function S ⊥ ( b, µ s , ν s ), and the beam functions , each evaluated at an appropriate scale to avoid large logarithms, and subsequently evolved to common scales. Here f j/N are the standard MS parton distribution functions and I gj and I gj are perturbatively calculable matching coefficients. Note that besides the usual renormalization scale µ, the soft function and beam function also depend on the rapidity scale ν. For p T m H the resummation is carried out by making the canonical scale choices which ensures there are no large logarithms in H, S ⊥ and B αβ g/Nγ in Eq. (3). Large logarithms are then resummed through the evolution factors U h (m H , µ B , µ h ) and U s (b, µ B , µ s ; ν B , ν s ), which connect the hard scale and the soft scales to the scales of the beam function, respectively. The U h and U s are derived from SCET renormalization group and rapidity renormalization group equations [39,40]. They depend on the cusp anomalous dimension Γ cusp [41], the gluon anomalous dimension γ V [42][43][44], the soft anomalous dimension γ s [45,46], and the rapidity anomalous dimension γ r [25,34], each of which is now known to three loops. In the case of the cusp anomalous dimension, the four-loop leading color approximation is also known [47]. The initial condition of evolution for the hard function and soft function are also known to three-loop order [42][43][44]46]. In the case of the beam function in Eq. (6), the initial conditions I gj and I gj are known to two loops [48,49], and the logarthmic terms are known to three loops, which involves the 3-loop splitting functions [41,50]. Our N k LL resummed calculation is obtained from these ingredients, where we include all logarithmic terms at O(α k s ) in the H, S ⊥ and B αβ g/Nγ boundary conditions. For larger p T ∼ m H the resumma-tion has to be turned off so that dσ r /dp 2 T = dσ s /dp 2 T in Eq. (1) and the cross section reduces to the fixed order result. This is achieved by making a transition of the various µ i to a single scale, and also transitioning to a single rapidity scale, using profile functions [51,52] as explained below.
The fixed-order (FO) distribution dσ f /dp 2 T is obtained as a perturbative expansion in α s by combining all parton-level contributions that contribute to a given order. As the transverse momentum is fixed to a nonvanishing value, the leading order process is given by the production of the Higgs boson recoiling against a parton [53]. At higher orders in perturbation theory, a method for the combination of infrared-singular real radiation and virtual loop corrections becomes necessary. In this letter, we employ the antenna subtraction method [54][55][56] to compute the fixed-order distribution up to NNLO in perturbation theory. The calculation is performed within the parton-level event generator NNLOJET [20], which combines the tree-level double real radiation corrections (RR, [57][58][59]), the one-loop single real radiation corrections (RV, [60][61][62]) and the two-loop virtual corrections (VV, [63]) with appropriate antenna subtraction terms. Distributions are obtained as binned histograms in p T . Schematically, the NNLO contribution to the distribution takes the form where the antenna subtraction terms dσ S,T,U N N LO (which upon integration add up to zero) ensure the finiteness of each n-parton phase space integral Φn . Details of the calculation are described in [20], where the trans-verse momentum distribution was computed to NNLO for large p T . Extending the lower bound on p T towards smaller values becomes increasingly challenging due to the large dynamical range probed in the phase-space integration and the associated numerical instabilities. We adapted the NNLOJET code to cope with this task and further split the integration region into several intervals in p T and applied dedicated reweighting factors in each region. With these optimizations, fixed-order predictions are obtained down to p T = 0.7 GeV, both for a linear binning of 1 GeV width and a logarithmic spacing of ten bins per e-fold.
At small transverse momentum p T with respect to the Higgs mass m H , the cross-section can be split into a singular (s) and non-singular (n) piece: with dσ s dp 2 dσ n dp 2 The coefficients c i,j are obtained analytically (up to the integrals over the PDFs) by setting the evolution factors U h and U s to unity in Eq. (3), and evaluating the hard function, soft function, and beam function at common scales µ i = µ F = µ R and ν s = ν B . To compare these singular terms with the full fixed order prediction, we integrate Eq. (10) over the same binnings that were used in the numerical evaluation in Eq.  Figure 1 compares the fixed-order contributions at LO, NLO and NNLO for the transverse momentum spectrum where the curve labeled as the SCET prediction is p T dσ s /dp T = 2p 2 T dσ s /dp 2 T at the corresponding order. For better visibility, the distributions are multiplied by p T , and each higher order contribution is displayed separately, instead of being added to the previous orders. The bottom panels show the difference of Comparison of the transverse momentum spectrum between fixed-order perturbation theory (FO) and singular terms from the expansion of the resummed prediction (SCET), using the sum of all partonic channels (pp) or with individual partonic channels (gg, qg, qq). In individual channels, q denotes the sum of quark and anti-quark of all flavors. the two curves, i.e. the non-singular parts which should behave as p T dσ n /dp T ∼ O(p 2 T ) for p T m H , to further elucidate the low-p T behaviour between the two predictions. The top frame shows that the small p T behaviour of the fixed-order spectrum is in excellent agreement with the predicted singular terms in Eq. (10). This agreement is further substantiated in the lower three NNLO Singular, |dσ s /dp T | NNLO Fixed Order, dσ f /dp T NNLO Non-Singular, dσ n /dp T

FIG. 2.
Comparison of full fixed-order spectrum, the absolute value of singular distribution, and the non-singular distribution through to NNLO. Here dσ n /dpT ∼ O(pT ) for pT mH .
frames, where individual parton-level initial states are compared (with q denoting the sum over quarks and antiquarks of all light flavours). We point out that the (numerically subdominant) qq channel turns out to be the numerically most challenging, since contributions from valence-valence scattering favor events with higher partonic center-of-mass energy than in any of the other channels. The excellent agreement between fixed-order perturbation theory and SCET-predictions for the singular terms serves as a very strong mutual cross check of both approaches. It demonstrates that our calculation of the non-singular terms is reliable over a broad range in p T , thereby enabling a consistent matching of the NNLO and N 3 LL predictions.
Matching and results.-For a reliable description of the transverse-momentum spectrum, the resummation of large logarithms in dσ s /dp 2 T has to be turned off at large p T . This can be seen clearly from Fig. 2, which depicts the full fixed-order spectrum, the absolute value of singular distribution, and the non-singular distribution, all through to NNLO. At p T 50 GeV, the singular distribution dominates the fixed-order cross section, and the resummation of higher order logarithms is necessary. Around 50 GeV, the singular and non-singular distribution become comparable, and resummation has to be gradually turned off. There are several different prescriptions on how to turn off the resummation [12,16,27,[66][67][68][69][70]. In this letter, we follow Ref. [16] by introducing b and p T dependent profile functions, defining where ρ(b, p T ) is used for µ s = µ s (b, p T ) = µ B , ν s = ν s (b, p T ), and µ h = µ h (p T ), which appear in Eq. (3). ρ l is the initial scale for each profile, taken to be the canonical scales in Eq. (7) so that at small p T the large logarithms are resummed. ρ r is the final scale for each profile, which is chosen to be µ h = µ B = µ s = µ F = µ R , while for ν s it is m H . The parameters s and t govern the rate of transition between the fixed order result and the resummation, where the transition starts at p T t − t/(2s), is centered at p T = t, and ends at p T t+t/(2s). In our calculation, we choose s = 1, and t = 20, 25,30,35,40,50 GeV to estimate the uncertainties from different profiles. The uncertainties for the final resummed + fixed-order prediction are estimated by three-point variations of i) the ρ l for µ h about m H and ρ r for all scales (varied simultaneously), and ii) the ρ l for µ B = µ s and ν s about b 0 /b (varied independently). We always fix ν B = m H . We take the envelope of the resulting 66 curves as the uncertainty band at each order. Further uncertainties in our calculation include the missing four-loop cusp anomalous dimension and the treatment of non-perturbative corrections at large b. They are estimated to be negligible compared with the aforementioned scale uncertainties. Additional independent uncertainties related to the parton distributions and value of α s (m Z ) should be included for a detailed phenomenological study. The final matched transverse momentum spectrum is shown in Fig. 3. We plot the distributions at LO+NLL, NLO+NNLL, and NNLO+N 3 LL. We also plot the unmatched NNLO distribution. At small transverse momentum, the fixed order distribution displays unphysical behavior, due to the presence of large logarithms. We see that the matched distribution smoothly merges into the fixed-order cross section around 40 GeV, and that the scale uncertainty reduces order-by-order in perturbation theory. The perturbative uncertainties at NNLO+N 3 LL have been reduced to < ∼ ±6% for 5 < p T < 35 GeV, are ±10% for intermediate p T , and decreasing again at large p T .
Conclusions.-In this letter we have presented for the first time precise predictions for the Higgs transverse momentum spectrum at small p T , with resummation at N 3 LL matched to fixed-order results at NNLO. The calculation builds upon efficient subtraction formalism for jet processes, improved formalism for resummation of large transverse logarithms, and known high-order anomalous dimensions and matching coefficients. We use an additive matching scheme which relies on the extraction of non-singular corrections from singular ones, and usually requires numerical Monte Carlo precision at the level of 1 per-mille, which imposes a strong challenge on fixed-order calculations in the infrared unstable small p T region. We have shown excellent agreement between SCET and NNLOJET in this region, which provides a highly nontrivial check of both calculations. The final matched predictions show a continuous reduction of scale uncertainties order by order, and are significantly more precise for small p T . We expect our results will have an important impact on understanding the detailed properties of the Higgs boson at the LHC. This calculation can be readily carried over to Drell-Yan production where high precision measurements are already available.
Acknowledgements.-We thank the University of Zurich S3IT and CSCS Lugano for providing the computational resources for this project. This research was supported in part by the UK Science and Technology Fa