Bottom-quark effects in Higgs production at intermediate transverse momentum

We provide a precise description of the Higgs boson transverse momentum distribution including top and bottom quark contributions, that is valid for transverse momenta in the range mb<pt<mt, where mb and mt are the bottom and top quark masses. This description is based on a combination of fixed next-to-leading order (NLO) results with next-to-next-to-leading logarithmic (NNLL) transverse momentum resummation. We show that ambiguities in the resummation procedure for the b-quark loops are of the same order as the related fixed-order uncertainties. We conclude that the current uncertainty in the top-bottom interference contribution to the Higgs transverse momentum spectrum is O(20%).

loops that facilitate the ggH couplings. For the processes gg → H + g, qg → Hq and qq → H + g these terms are sensitive to gluons emitted from both "inside" and "outside" the loops at finite transverse momentum, i.e. to the structure-dependent radiation.
Light-quark contributions to Higgs production in gluon fusion make the resummation of the transverse-momentum distribution difficult [31,32]. Indeed, since both top and bottom quark loops contribute to the ggH coupling and since these loops are characterized by very different intrinsic scales for the structure-dependent radiation (m b and m t ), it appears that one will have to treat them differently. However, this is not possible since the dominant contribution is given by the interference of the two amplitudes. In addition, since it is not understood how to resum the potentially large logarithms log(p ⊥ /m b ) that appear in the light-quark loops, it becomes impossible to treat all the different contributions to the Higgs p ⊥ spectrum on the same footing. The best thing that one can do is to employ a variety of prescriptions for combining light quark contributions with smallp ⊥ resummations and to study how the resulting uncertainty in predictions compares with other sources of theoretical error.
The goal of this paper is to study the Higgs p ⊥ spectrum including top and bottom-quark contributions at next-to-leading order combined with next-to-next-to-leading logarithmic transverse momentum resummation (NLO+NNLL). A similar study at leading order combined with next-toleading logarithmic resummation (LO+NLL) was performed in Ref. [31]. 2 To this end, we include the recently computed NLO QCD corrections to light-quark contributions to Higgs production in gluon fusion [29,30,33]. We find that the uncertainty in our matched NLO+NNLL result for the top-bottom interference contribution to the Higgs transverse momentum distribution in the region 10 GeV < ∼ p ⊥ < ∼ 100 GeV is dominated by ambiguities in the perturbative description of light-quark loops rather than by uncertainties in the resummation itself. In particular, we do not find large uncertainties related to the choice of the resummation scale for the b-quark loops.
The paper is organized as follows. In Section 2 we briefly review the structure of small-p ⊥ resummation for the case of point-like interactions, and elucidate its main assumptions and limitations. We also study light-quark contributions, discuss why in this case the resummation is challenging and describe a possible pragmatic solution to this problem. In Section 3, we explain the implementation of the resummation procedure for the b-quark contribution and study its ambiguities, and we present our main results for the Higgs transverse momentum distribution. We conclude in Section 4. Some useful formulas and derivations are collected in the Appendix.

The standard point-like case
We would like to describe the transverse momentum distribution of Higgs bosons produced in hadron collisions. This is non trivial and requires a combination of fixed order and resummed perturbative calculations. Indeed, depending on the value of the Higgs boson transverse momentum, we can distinguish two regions. For large values of transverse momenta p ⊥ ∼ m H , one can compute dσ/dp ⊥ in a perturbative expansion in α s following standard rules of perturbative Quantum Field Theory. For small values of the transverse momentum p ⊥ m H , the situation is different since emerging large logarithms ln(p ⊥ /m H ) 1 may compensate the smallness of the strong coupling constant, α s ln 2 (p ⊥ /m H ) ∼ 1, and spoil a conventional perturbative expansion. To deal with this case, one resums the logarithmically enhanced terms to all orders in the coupling constant, and develops a perturbative expansion on top of the resummed result.
Since, eventually, we need to describe the Higgs boson p ⊥ distribution for all values of transverse momenta, the two distinct approaches -resummation and fixed order computations -have to be combined. This is done by smoothly interpolating between results derived at small and large p ⊥ . The region where the transition happens is characterized by a quantity that we refer to as the resummation scale Q. This scale has the following physical meaning: for p ⊥ < ∼ Q, the transverse momentum distribution is mostly described by the resummed result, while for p ⊥ > ∼ Q it is mostly described by the fixed order computation.
In order to discuss these concepts more precisely, we consider the all-order resummation in a toy model, where we work at leading-logarithmic (LL) accuracy. To this end, we consider the cumulative distribution Σ(p ⊥ ) = p ⊥ 0 dp ⊥ dσ dp ⊥ . (2.1) At low p ⊥ , we resum the logarithms of ln p ⊥ /m H and write In this region, the distribution is dominated by the emission of soft and collinear partons. In the LL approximation it is sufficient to consider the most singular contribution to the QCD matrix elements, where all final-state partons are soft and strongly ordered in angle. In this limit, the squared matrix element for the emission of n extra partons gg → H + n is given by the product of the matrix element for gg → H times n independent eikonal factors. More specifically, at LL the partonic p ⊥ distribution can be simplified as dσ dp ⊥ [dp H ]|M(p 1 + p 2 → H)| 2 δ (4) (p 1 + p 2 − p H ) where σ 0 is the Born cross section for gg → H. The overall exponential factor contains the all-order effects of soft-collinear virtual gluons which are encoded in the leading divergence of the gluon form factor M(p 1 + p 2 → H). 4 The distribution in the small p ⊥ region is governed by two competing mechanisms. In the strict limit p ⊥ → 0, the dominant contribution comes from emissions with finite transverse momentum p ⊥ k ⊥i m H that mutually cancel in the transverse plane. This collective effect gives rise to a power suppressed scaling [34] As p ⊥ increases, but still remains small compared to m H , the distribution is described by kinematic configurations with p ⊥ ∼ k ⊥i m H . As discussed in Appendix A.1, in this region the cumulative distribution features an exponential suppression of the form At larger transverse momenta (p ⊥ ∼ m H ) the approximation that led to Eq. (2.5) is not justified anymore. Therefore, in this region one has to smoothly switch from the resummed prediction to the fixed-order one, where the effect of the hard radiation is treated correctly. This can be done for example using the following matching formula where we indicate with T f.o.
[f ] the fixed-order expansion of f . At small p ⊥ the difference between the fixed-order result and the Taylor expansion of the resummed result is free of logarithmicallyenhanced terms lim This allows one to extend the fixed-order description to p ⊥ → 0 and, at the same time, ensures that all terms that contain large logarithms at low p ⊥ are resummed. The precise way to switch from the resummation to the fixed-order description is ambiguous. One source of ambiguity comes from choosing a particular form for the matched cross section (in our example, Eq. (2.8), we chose to combine the resummed and fixed-order predictions additively). A second source of ambiguity is connected with the scale at which the transition from resummed to fixed-order result takes place. Although all of these effects are formally of higher-order both in the resummation and fixed-order counting, their numerical impact can be non-negligible. We consider the latter issue in what follows, while leaving a discussion of the choice of the matching scheme to the next section.
In order to switch off resummation effects at large p ⊥ , one can modify the resummed cross section by including controlled power-suppressed corrections. One possible way to do this is to modify the resummed logarithms in Eq. (2.7) as follows 5 where Q is an arbitrary scale of order m H . Moreover, we write where p is a positive number. The motivation for the transformations described above is as follows: • First, we split the resummed logarithm L into the sum of a small logarithm ln(m H /Q) (with Q ∼ m H ) and a large logarithm ln(Q/p ⊥ ). This operation allows us to introduce a generic scale Q which then appears in the resummed result. We can now expand L around ln(Q/p ⊥ ), retaining all terms with the desired logarithmic accuracy. Effectively, this implies that ln(m H /Q) is treated perturbatively at fixed order. In our LL example, for p ⊥ ∼ k ⊥i m H , this means where all terms beyond LL were neglected. This prescription is convenient because the Qdependence is always of higher-logarithmic order and, therefore, a Q-variation probes the size of subleading logarithms that are not considered in the resummation.
• Second, we modify the logarithm ln(Q/p ⊥ ) by including power-suppressed terms that forceL to vanish at large p ⊥ . These modifications do not affect the small-p ⊥ limit. Indeed, it follows from Eq. (2.11) that As a consequence, the resummation scale Q and the scaling parameter p must be chosen in such a way that the high-p ⊥ scaling of the resummed component (and its fixed-order expansion) does not modify the scaling of the fixed-order prediction. This means that p and Q are to be chosen in such a way that the resummed component vanishes more quickly than the fixed-order result for p ⊥ > ∼ Q.
The above discussion shows that Q is indeed the scale at which the transition between resummed and fixed-order results occurs. Similarly to the renormalization and factorization scales, its choice is ambiguous, although certain conditions should be satisfied. Indeed, it is clear that (a) Q should not be too different from m H , to ensure that ln Q/m H are not large and (b) it should approximately correspond to the scale at which the soft and collinear approximations to the matrix element and kinematics break down. In practice, one can choose Q by comparing the exact result Σ f.o. with the expansion of the resummed result T f.o. [Σ resum ], and set Q to the p ⊥ scale at which the two start to significantly deviate from each other. This is illustrated in Figure 1, which shows the difference between the LO differential p ⊥ spectrum and the expansion of the resummed result at the same order. Specifically, we plot (2.14) We observe that when only the top contribution is included (solid, red curve), the logarithmic terms account for about half of the fixed-order result at scales p ⊥ ∼ 50 − 60 GeV. This suggests that the resummation scale should be of this order. We conventionally choose Q = m H /2 as a central value. As far as the choice of the parameter p is concerned, we have to ensure that at large p ⊥ the resummed component vanishes faster than the fixed order. Considering the asymptotic scaling in Eq. (2.13), we choose p = 4 which guarantees that the differential distribution vanishes as fast as 1/p 5 ⊥ for p ⊥ > ∼ Q. In principle, any value of p greater than 3 will equally do, since p only determines how fast the resummation is turned off above the scale Q. We have indeed checked that by varying p by one unit around p = 4 the results do not change significantly. /d p top only top and bottom interference The same figure shows results for the full spectrum where both top and bottom loops are included (dotted, blue curve), and results for top-bottom interference (dot-dashed, green curve). 6 In this case the situation changes considerably, and this will be the subject of the next section.

Issues with b quarks
The "standard" approach to resummation described in Section 2.1 becomes problematic in case of the Higgs boson production in hadron collisions [31,32,35]. The difficulty is related to the fact that the ggH vertex is not point-like but, rather, is induced by a quark loop. The presence of the quark loop implies the existence of structure-dependent radiation with peculiar properties and has important consequences for the resummation. The key to the following discussion is the appreciation of the fact that the structure-dependent radiation is suppressed if p ⊥ is smaller than the mass of the quark but it becomes important otherwise. For p ⊥ larger than the quark mass, the soft and collinear approximations that provide the foundation for small p ⊥ resummation become unreliable, as they focus on emissions off external lines and systematically neglect structure-dependent effects. In this section we elaborate on this issue.
We consider Higgs boson production in gluon fusion mediated by a quark loop. We denote the mass of the quark by m q and consider two cases m q m H and m q m H . In the first case, the structure dependence enters at p ⊥ > ∼ m q m H , so that emissions off the external lines dominate for transverse momenta up to the Higgs mass and even higher. Therefore, if we restrict ourselves to values of p ⊥ that are comparable to m H , the situation is no different from point-like interaction, and there are no issues in the resummation procedure described in the previous section. In the Standard Model, this is indeed what happens with the top loop contribution to the Higgs boson transverse momentum spectrum.
The second case, m q m H is very different. Indeed, in this case there are three distinct regions In the first region p ⊥ < ∼ m q , the transverse momenta of the Higgs boson and the recoiling partons are typically small enough not to resolve the structure of the loop and the extra radiation factorizes. For the bottom-quark contribution (m q ∼ 5 GeV) the effect of additional QCD radiation is strongly suppressed in this region by all-order effects, so that its impact on the total cross section is small. 7 Note, however, that in this region there are large logarithmic contributions of the type ln 2 m H /m b , whose resummation is not fully understood even at the lowest perturbative order [36,37].
In the second region m q < ∼ p ⊥ < ∼ m H the structure-dependent radiation becomes essential and the ggH vertex does not factorize. In addition to the usual logarithms ln m H /p ⊥ , the radiation gives rise to logarithms ln p ⊥ /m q and ln m H /m q , whose origin and potential resummation are not well understood. 8 The reason why the small-p ⊥ resummation is problematic in this region is the following. Emissions off internal lines can become as important as emissions off external lines, and together they probe the loop structure of the ggH vertex. It follows that approximating the small p ⊥ region with an on-shell ggH form factor is not justified. In particular, while form factor effects in the top-quark case only introduce (p ⊥ /m H )-suppressed corrections, in this case they both introduce a new logarithmic structure (ln p ⊥ /m q and ln m H /m q ) and suppress radiation with p ⊥ > ∼ m q . In other words, while in the top case, described in Section 2.1, at finite p ⊥ , the coefficients of the logarithms differ from the resummed result by p ⊥ -suppressed terms, in the b-quark case, this difference contains new logarithmic terms ln p ⊥ /m b and ln m H /m b in the region m b < ∼ p ⊥ < ∼ m H . As a consequence, the collinear approximation should not be expected to work far away from the b-quark threshold. To quantify this effect, we go back to Figure 1. We see that, while for the top-only case (solid, red line) the collinear approximation to the leading order accounts for half of the result at about p ⊥ ∼ 50 − 60 GeV, for the top-bottom interference this scale is reduced to about 30 GeV. When top and bottom contributions are considered together, this effect is less dramatic since in the SM the interference accounts for about ∼ 5% of the full result. This can be seen in the dotted blue curve of Figure 1.
Because of the above issues, it is clear that constructing a reasonable description of the b-quark contribution to the Higgs transverse momentum distribution is problematic. Since, as we already stressed, in this case the resummation of potentially large logarithms is not entirely understood, the best we can do is to use different ways to interpolate between regions of small and large transverse momenta and check to what extent the different results are compatible.
As already stated, the Higgs boson production in the Standard Model is dominated by the top quark loop; the bottom loop provides a very small contribution that is lifted up to O(−5%) by its interference with the top amplitude. Because of this, a O(20 − 30%) control on the top-bottom interference is sufficient to control the Higgs transverse momentum spectrum at the few percent level. With this in mind, we now study in more detail the different ways to treat the bottom contribution.
One option is to apply Eq. (2.8) with the resummation scale set to Q ∼ m b [31]. This choice is equivalent to employing fixed-order description for all values of transverse momenta. Indeed, for 1 GeV < ∼ p ⊥ < ∼ Q ∼ m b , the resummed logarithms ln(Q/p ⊥ ) never become large and for p ⊥ > ∼ Q the fixed-order result is adopted anyhow. Since a typical error made within this approach is provided by uncalculated higher order terms, if we use a NLO computation for the interference, we make an error of order 9 The previous option amounts to neglecting the resummation for the top-bottom interference and to using the fixed-order result for all transverse momenta; the other extreme alternative consists of extending the resummation beyond its established domain of validity. We can do this by using the same resummation scale, Q ∼ m H /2, both for the top and the top-bottom interference contributions [32,39]. In this case, at low p ⊥ we introduce logarithms ln p ⊥ /Q in the interference through the resummation prescription which are not guaranteed to be correct and, by doing that, we effectively introduce errors that are similar to those discussed above. At higher p ⊥ the impact of these logarithms becomes smaller since at p ⊥ ∼ Q the resummation effects smoothly turn off and we recover the fixed order prediction. Hence, when we under-or over-resum logarithms we expect comparable O(20%) theoretical errors on the interference contribution to the Higgs p ⊥ spectrum. An important question is whether these different sources of uncertainties pull the predictions apart or they remain compatible with each other.
Before concluding this section, we mention that within the additive matching scheme of Eq. (2.8) the resummation term, which is proportional to the lowest-order form factor at zero transverse momentum, is added to the fixed order result. As we mentioned earlier, the form factor effects lead to a dependence of the leading order amplitude on p ⊥ , that is not captured in this approach. To account for this, we also consider a multiplicative matching scheme, which can be schematically defined as 10 Similarly to Eq. (2.8), Eq. (2.15) smoothly interpolates between a low p ⊥ Q region, where resummation dominates, to a large p ⊥ Q region, where the result is obtained from a fixed order calculation. Clearly, the fixed order accuracy is preserved in the p ⊥ → 0 limit. The main difference with Eq. (2.8) is that now the higher order terms induced by the resummation in the transition region are weighted with the fixed (lower) order result at finite transverse momentum. This should at least partially capture the p ⊥ dependence of the exact higher order amplitude, and lead to a more realistic description of the physics. Because of this, we choose the multiplicative matching scheme Eq. (2.15) as our default matching scheme. Nevertheless, matching ambiguities are by construction of higher-order nature and, therefore, any matching prescription is formally equally valid. Differences between matching prescriptions can be used to estimate the uncertainty in the transition region.

Inclusion of bottom-quark loops and matching uncertainties
In this section we describe the practical and technical implementation of the top-bottom interference in the resummation and matching, and the uncertainty associated with it. As we described in Section 2.2, the rigorous resummation in the presence of the bottom-quark loop is currently impossible. To remedy this problem, we adopt different approaches to include this contribution in the matched result. We use the arbitrariness in the choice of the resummation scale associated with the top-bottom interference and in the choice of the matching scheme to assess the inherent ambiguity of the resummation procedures.
We start by discussing the resummation scale. We treat separately the contribution of the top-squared amplitude and the top-bottom interference. 11 In particular, we associate two different resummation scales with the top and the interference contributions, and we use the following notation to denote the various cumulative distributions As explained in Section 2.1, for the top-only contribution we set the resummation scale to Q t = m H /2. For the interference, instead, we use the following prescriptions to quantify the associated uncertainty (see Section 2.2): • We choose Q b ∼ m b , effectively switching off the resummation for the interference at scales of the order of the bottom mass. As it was initially suggested in Ref. [31], we choose Q b = 2m b as our central scale. This is achieved by computing This implies that in the region of transverse momenta that we are interested in, the interference is described only at fixed order and no resummation for this contribution is performed.
• We consider the opposite situation in which we rely on the collinear approximation also for m b p ⊥ , and simply treat the new logarithmic terms that appear above this scale as a regular remainder that can be described at fixed order. As a consequence, the resummation for the interference contribution is switched off, as in the top-only case, at scales of order 60 GeV. We choose Q t = Q b = m H /2 as our central scale, for simplicity.
In both approaches, logarithms of the ratio p ⊥ /m b are not resummed. Although in the region m b p ⊥ m H these logarithmic terms can be potentially large and therefore should be included to all orders, recent studies seem to suggest that an accurate prediction of these terms is achieved by considering the first few terms in the fixed-order perturbative expansion [36,37].
As far as the resummation is concerned, the result will be nearly identical to the m q → ∞ one. 12 The only difference is that now the Born squared amplitude and the hard-virtual correction will contain the full dependence on the top and bottom masses. In particular, no modification of the p ⊥ -dependent radiation pattern is introduced. Technically, we implement the LO and the NLO amplitudes for gg → H with full mass dependence following Ref. [40].
We now study numerically the difference between the two prescriptions for the bottom resummation scale. We start by introducing the setup that we adopt for our predictions. We consider proton collisions at the 13 TeV LHC. The Higgs boson mass is taken to be m H = 125 GeV and the top and bottom pole masses 13 are set to m t = 173.2 GeV and m b = 4.75 GeV, respectively. We work within a fixed flavor-number scheme (n F = 5) and use the PDF4LHC15_nnlo set of parton distribution functions [41] interfaced through LHAPDF6 [42]. We use the value of the strong coupling constant α s provided by the PDF set. As central values for the renormalization and factorization scales we take In order to estimate the perturbative uncertainties for each prediction, we perform a 7-point variation of the factorization (µ F ) and renormalization (µ R ) scales around the central value by a factor of two. Moreover, we vary Q t and Q b by a factor of two around their respective central values, keeping fixed µ R = µ F = M T /2. The final uncertainty band is obtained as the envelope of all above variations. As a default, we adopt the multiplicative scheme discussed in Section 2.2 and described in detail in Appendix A.3.
The fixed-order NLO results for the top-bottom interference are based on the calculation presented in [33], which in turn comprises the two-loop amplitudes for gg → Hg, qq → Hg and qg → Hq derived in [29,30] together with corresponding loop-squared real radiation amplitudes as provided by OpenLoops [43,44] combined with Collier [45]. For the Monte Carlo integration and subtraction the Powheg-Box-Res is used [46,47]. 12 HEFT in the following. 13 We work in the on-shell renormalization scheme as a default.
We now discuss the dependence on the choice of the resummation scale associated with the bottom contribution. We start by comparing results for the top-bottom interference for two values of the resummation scale Q b . The results are displayed in the left plot in Figure 2. The two -0.14 -0.12 On-shell, multipl. scheme  predictions differ by about 20% for p ⊥ ∼ 30 GeV, in line with what we expected from the discussion in Section 2.2. We note, however, that although the two results are computed for very different choices of the resummation scales, they are still compatible within their respective uncertainties. Only for p ⊥ > ∼ 50 GeV the two results differ significantly. However, in this region the contribution of the interference is already very small. Each of the two results has a relative uncertainty of about 15% for p ⊥ < ∼ 40 GeV. The variations of the resummation scales around their central value, and the variation of µ R and µ F have a similar impact on the final band.
The right plot of Figure 2 shows an analogous comparison for the transverse momentum distribution that includes both top and bottom contributions. Since the interference only accounts for about 5% of the full result, we find that the two resummation prescriptions for the top-bottom interference are indistinguishable within the uncertainties of the top contribution. Indeed, in this case the uncertainty band is dominated by the µ R and µ F variation of the top contribution, which amounts to about 10 − 15% for p ⊥ < ∼ 40 GeV, while the resummation-scale uncertainty amounts to about 5% in this region. Note however that the top-only contribution has been computed one order higher, both in fixed-order QCD [7][8][9][10][11] and in the resummation framework [26]. In this paper, we focus on the b-quark effects and hence do not include these results but, as a matter of principle, they can be used to further reduce the uncertainty on the top contribution.
We now investigate the second source of resummation ambiguity, namely the choice of the matching scheme. As discussed in Section 2.2, besides our default multiplicative scheme we also consider an additive scheme. Both schemes are precisely defined in Appendix A.3. We remind the reader that, as far as the top-bottom interference is concerned, the main qualitative difference between the two approaches is that within the additive matching scheme, the resummation contribution is proportional to the gg → H form factor at zero transverse momentum whereas in the multiplicative matching scheme it is weighted by the form factor g * g * → H at finite transverse momentum. In order to study this source of ambiguity more precisely, we consider the additive matching scheme, with two different scales (Q b = m H /2 and Q b = 2m b ), and compare the results to the multiplicative scheme.
Since in the additive scheme the resummed contribution does not include form-factor effects,  we expect sizable differences between results obtained with Q b = m H /2 and Q b = 2m b . We recall that this is not the case in the multiplicative scheme (see Figure 2) where form-factors effects are automatically accounted for. We then show, in Figure 3, the comparison between the top-bottom interference in the default multiplicative scheme with Q b = m H /2 and the additive scheme with Q b = 2m b (left plot) and Q b = m H /2 (right plot). We observe that the difference between the two schemes is larger when the additive scheme with Q b = m H /2 is used. Nevertheless, we find that also in this case, the difference between the two schemes for the interference does not exceed ∼ 20% in the bulk of the distribution. Again, the full transverse momentum distribution, shown in the left plot of Figure 4, is only mildly affected by this ambiguity. Finally, in the right plot of Figure 4, we show the ratio of the full distribution computed using the default multiplicative scheme, to the corresponding HEFT result. The default result, i.e. multiplicative matching scheme with Q t = Q b = m H /2, is in good agreement with the NLO prediction. For comparison, we also report the other extreme solution obtained with the multiplicative scheme with Q t = m H /2, Q b = 2m b . We observe that this choice is in good agreement with both the fixed order and the default matched solution.
In summary, we find that a conservative approach towards the inclusion of bottom-mass effects in the matched prediction for the Higgs-p ⊥ spectrum leads to a ∼ 20% ambiguity on the top-bottom interference in the region m b < ∼ p ⊥ . Since the interference provides a rather small contribution to the Higgs transverse momentum distribution, this ambiguity translates into a few-percent uncertainty in the Higgs p ⊥ spectrum.
In what follows, we will use the result obtained with the multiplicative matching scheme with Q t = Q b = m H /2 as our central value. To estimate uncertainties, we will consider the envelope of scale variations and the result obtained either by using the multiplicative or the additive scheme with Q t = m H /2, Q b = 2m b . In addition to these source of uncertainty, an additional ambiguity arises from the choice of the renormalization scheme for the quark masses. This will be discussed in the next section.

Mass-scheme uncertainty and final results
In this section, we present our final results for the NNLL+NLO matched distributions. We use as default the multiplicative matching scheme with resummation scales Q t = Q b = m H /2. We renormalize the bottom-quark mass in the on-shell scheme. To estimate the uncertainty we change the details of the resummation and matching as explained in the previous section. In addition,  we consider a different renormalization scheme for the bottom quark mass to estimate the related uncertainty. To this end, we employ the MS renormalization scheme. We take the mass renormalization scale to be 100 GeV, and use m b = m MS b (100 GeV) = 3.07 GeV as an input parameter. 14 In Fig. 5 we display the results for the top-bottom interference contribution. The fixed order result is presented in the left plot. We show the uncertainty band for the on-shell mass-renormalization scheme and the central value for the MS scheme. The uncertainty band is calculated from a 7-point scale variation. We see that the scheme ambiguity is larger than the scale variation, as already observed in Ref. [33]. The right plot shows our results for the matched distributions with the two different mass schemes. The difference between the two bottom-mass schemes is similar to the fixed order case, but since now the matched prediction has a smaller uncertainty, the separation between the two results is more significant. As follows from  has an ambiguity of about 15 − 20% down to p ⊥ ∼ 10 GeV. In order to improve on this, one would need a NNLO calculation for the top-bottom interference, which is currently out of reach.
The analogous plots for the full distribution that includes both top and bottom amplitudes are shown in Fig. 6, for the fixed order (left) and resummed (right) results. Unlike for the top-bottom interference contribution, in this case the difference between the two results for the bottom-mass schemes are much smaller, at the level of a few percent. This is because the top-bottom interference contributes to just O(−5%) of the full spectrum.
Our best current predictions for the top-bottom interference and full p ⊥ spectrum including all the relevant uncertainties are shown in Fig. 7. As discussed earlier, the uncertainty bands are obtained as an envelope of: • a 7-point renormalization and factorization scale variation; In addition, if the default matched result but in the MS renormalization scheme for the bottomquark mass is outside the uncertainty band estimated as described above, we extend the uncertainty band to accommodate the mass scheme ambiguity. In fact, as shown in Figure 5, the latter ambiguity is the major source of uncertainty for the top-bottom interference for transverse momenta below 30 GeV.
The top-bottom interference is shown in the left plot of Fig. 7. The qualitative features of the fixed-order result are unchanged by the resummation, which however has a noticeable effect on the shape of the distribution. Our final result has an uncertainty of about ∼ 20%, and is compatible with the fixed-order one. In the right plot of Fig. 7 we present the results for the full spectrum. At large values of the Higgs p ⊥ > ∼ 30 GeV the fixed order result is contained in the error band of the resummed result. However at smaller values, p ⊥ < ∼ 30 GeV we observe a marked difference between the two results. The error for the full matched result is close to 10% for p ⊥ < ∼ 30 GeV and close to ∼ 20% at larger p ⊥ . We stress however that the uncertainty on the dominant top contribution can be further reduced by employing the results of Refs. [7-11, 26, 48].

Conclusions
In this paper we performed a detailed study of the Higgs transverse momentum distribution, focusing on the region of intermediate values of transverse momenta, Indeed, a precise theoretical control of the Higgs p ⊥ distribution in this region is essential to test the Higgs sector of the Standard Model. In particular, it provides a rare opportunity to probe the Yukawa couplings of light quarks, which are currently poorly constrained. In fact, although the main contribution to the Higgs production cross section is due to the coupling of the Higgs to top quarks, the coupling to bottom quarks has a non-negligible impact on the total cross section through its interference with the top, decreasing the cross section by about O(5%).
The theoretical description of the Higgs p ⊥ distribution for m b < ∼ p ⊥ < ∼ m H in QCD is particularly challenging since, once the contribution of bottom quarks is included, the perturbative cross section for small p ⊥ suffers from the presence of potentially large logarithms ln (p ⊥ /m b ), ln (m H /m b ), which can spoil the convergence of the perturbative expansion. The physical origin of these large logarithms is not yet fully understood, and their all-order resummation remains currently out of reach.
Given these conceptual limitations, we provided our best theoretical description of the Higgs p ⊥ distribution at NNLL+NLO QCD for moderate values of the transverse momentum, including dependence on the bottom mass. An important part of our study was a proper assessment of the theory uncertainty of our results. The NLO result for the top-bottom interference suffers from scale uncertainties, which amount to around 15%. On top of this, a non-negligible source of uncertainty is provided by the renormalization scheme ambiguity for the bottom-quark mass, which we estimated by varying from the on-shell to the MS scheme. This amounts to an uncertainty of up to 20% and it dominates the error budget of our prediction for the top-bottom interference at small values of the Higgs p ⊥ . Together with the uncertainties associated with the fixed order calculation, we also performed a detailed study of the ones associated with the resummation procedure in the presence of bottom quarks. In order to estimate these ambiguities for the top-bottom interference, we matched the fixed order NLO predictions with the NNLL resummed cross-section using two different schemes, an additive and a multiplicative one, and two very different choices of the resummation scale, Q b = 2m b and Q b = m H /2. This leads to an uncertainty between 15−20% on the top-bottom interference contribution to the p ⊥ spectrum. Since the interference amounts to about 5% of the full p ⊥ spectrum, we conclude that unknown higher order b-quark mass effects can modify the Higgs transverse momentum distribution by few percent. All ambiguities associated with the resummation in the presence of bottom quarks produce consistent results within the NNLL+NLO uncertainty band, which is however driven by uncertainties in the (NLO) top quark contribution. The latter is currently known to higher N 3 LL+NNLO accuracy [7-11, 26, 48]. It would be interesting to combine these results with the ones presented in this article. We leave this for future investigations.
In conclusion, we presented a description of the Higgs p ⊥ spectrum at NNLL+NLO QCD including both top and bottom quark contributions. We found that the uncertainty on the top-bottom interference is O(20%) in the region of interest m b < ∼ p ⊥ < ∼ m H . Given the intrinsic ambiguities from scale dependence and, in particular, from the choice of the bottom-mass renormalization scheme and matching scheme, any improvement in this description will inevitably require the computation of the NNLO QCD corrections to the bottom-quark contribution to gg → H and gg → H + jet.
Moreover, if k 1 is collinear top 1 one has (from (1 − z An analogous limit on z where the product runs over all emissions off leg 1 that occur prior to the emission we are parametrizing, and we used the fact that in the soft limit z A similar parametrization holds for emissions off leg 2 .
The transverse recoil of the radiation is absorbed entirely by the Higgs boson that acquires a transverse momentum In order to predict the p ⊥ → 0 limit, we need to sum emissions at all orders in the strong coupling.
With LL accuracy, the squared amplitude for n emissions can be approximated by a product of n independent splitting kernels, as the soft correlation between emissions starts contributing at NLL order. The physical picture corresponding to this approximation is given by a set of independent emissions off legs 1 and 2 . In this approximation, the differential partonic distribution can be written as where the eikonal squared amplitude for a single emission reads In Eq. (A.7), the coupling is evaluated at k ⊥ to account for the leading-logarithmic contribution of the gluon branching into either a pair of soft quarks or gluons, see e.g. [51] for a detailed explanation. The resummation is naturally performed at the level of the cumulative distribution, defined as Indeed while the differential spectrum involves plus distributions in p ⊥ , Σ(p ⊥ ) is a regular function. From Eq. (A.6), it follows that the cumulative distribution with LL accuracy can be written as where f g (µ F ) is the gluon parton density evaluated at the factorization scale µ F , and the convolution is defined as usual Since p ⊥ only constrains the transverse momentum of the emissions, we can perform the integrals over the z ( i) i components inclusively. It is therefore convenient to introduce the functions ) . (A.11) This notation allows us to parametrize the real-emission matrix element and phase space as where we defined ζ i = k ⊥i /k ⊥1 . We now discuss the purely virtual corrections, which are encoded in the gluon form factor |M(p 1 +p 2 → H)| 2 . We write it as where the function H contains all the IRC singularities and the constant finite corrections of the form factor, and M B denotes the Born amplitude. Since we are working with LL accuracy, we are only interested in the leading singular term of H at all orders (while neglecting all finite terms) which can be written as Note that the integral in Eq. (A.14) is divergent and is to be considered as regularized. In order to cancel the IRC divergences of the real emissions (A.6) against the ones in the virtual corrections (A.14) at all orders, we introduce a small slicing parameter > 0 such that all emissions with a transverse momentum k ⊥i smaller than k ⊥1 can be ignored in the computation of the observable p ⊥ , in the limit → 0. The real emissions with k ⊥i < k ⊥1 , hereby denoted as unresolved, can be directly combined with the virtual corrections at all orders. Their combination gives rise to an exponential suppression factor of the type On the other end, emissions with k ⊥i > k ⊥1 , that we denote as resolved, are constrained by the observable's measurement function and therefore cannot be integrated over inclusively. The resummed LL cross section thus reads where σ B is the Born cross section. The above formula, in the limit → 0 exactly reproduces the LL corrections to the p ⊥ distribution, see Ref. [26] for a formal proof. Eq (A.16) can be further simplified by observing that in the resolved radiation one always has ζ i ∼ 1, since configurations in which ζ i 1 are automatically canceled against the exponential Sudakov factor e −R( k ⊥1 ) . Therefore, one can expand the functions R (ζ i k ⊥1 ) in powers of ln(1/ζ i ) as .17) and retain terms that contribute at a given logarithmic order. In particular, at LL, only the first term in this expansion contributes, and higher-order terms matter at higher logarithmic orders (see Refs. [26] for details). Similarly, we can consistently expand out the dependence of the exponential Sudakov as where the dependence manifestly cancels against the one in the resolved contribution, and we defined Therefore, with LL accuracy, Eq. (A.16) becomes Equation (A.20) is suitable for a numerical implementation, as explained in Ref. [26] in detail. The dependence on is at most power suppressed (i.e. O( p ⊥ )) and it vanishes in the limit → 0. This limit can therefore be taken safely numerically, and the result is absolutely stable for very small values of . 16 We now introduce the resummation scale Q as a possible way to switch off the resummation at large transverse momentum. This is defined with a procedure similar to the one discussed in the text. We first break the logarithm as follows The above operation will allow us (as explained shortly) to have an additional handle (namely the scale Q) to estimate the size of subleading logarithmic terms. Moreover, we also slightly modify the phase space available for the radiation, by introducing power-suppressed contributions that ensure that at large p ⊥ the resummation effects completely vanish. This can be done, as a first step, by modifying the resummed logarithms as follows where p is a positive real parameter which is chosen such that the resummed differential distribution vanishes as 1/p p+1 ⊥ at large p ⊥ . The above prescription essentially amounts to the following 1. First, we split the resummed logarithm L into the sum of a small logarithm ln(m H /Q) (with Q ∼ m H ) and a large one ln(Q/k ⊥1 ). This operation allows one to introduce a generic scale Q which appears in the resummed logarithms. One can now expand L about ln(Q/k ⊥1 ), retaining all terms with the desired logarithmic accuracy. Effectively, this implies that ln(m H /Q) is treated perturbatively at fixed order. Moreover, we replace ln(Q/k ⊥1 ) by the modified logarithmL. In our LL example this means whereR andR are functions of the modified logarithmL only.
This corresponds to the Jacobian for the transformation (A. 22), and ensures the absence of fractional (although power suppressed) α s powers in the final distribution [26]. This factor, once again, leaves the small k ⊥1 region untouched, and only modifies the large p ⊥ region by power-suppressed effects. This is effectively mapping the limit k ⊥1 → Q onto k ⊥1 → ∞.
Although this procedure seems a simple change of variables, we stress that the observable's measurement function (i.e. the Θ function in Eq. (A.20)) is not affected by this prescription. As a consequence, the final result will depend on the parameter p through power-suppressed terms.
The difference between the above prescription and what was introduced in the text is that the argument of the (modified) logarithms is now k ⊥1 instead of p ⊥ . This prescription is technically more correct, since in the small k ⊥1 region, which governs the p ⊥ → 0 limit, the modified logarithms leave Eq. (A.20) untouched. Conversely, at large k ⊥1 , where one has k ⊥1 ∼ p ⊥ , the above prescription reduces to what was defined in the text, i.e. the modified logarithms of k ⊥1 in this region are formally equivalent to modified logarithms in p ⊥ . To see this, we observe that when k ⊥1 Q the function R (k ⊥1 ) 1. Therefore, the probability of having any emission after the first one in Eq. (A.20) is strongly suppressed. As a consequence, at large k ⊥1 , the only relevant event is the one that involves a single emission k 1 , for which the cross section reads It is easy to see that, if Eq. (A.25) were evaluated without the factor J , it would lead to additional power-suppressed terms with fractional power of the coupling, which are clearly spurious.

A.2 Final formulas for NNLL resummation
Beyond LL, Eq. (A.20) is corrected to account for the description of the real-emission matrix element and phase space in less singular configurations, as well as higher perturbative corrections. To NNLL order it can be expressed as [26] where ζ i = k ⊥i /k ⊥1 . In this formula, the Sudakov radiatorR(k ⊥1 ) is corrected with respect to its LL expression by higher-order corrections of both soft and collinear origin. The same comment applies to the functionR i which we decided to split into the oldR (derivative of the LL radiator defined above), plus a correction that contains all subleading effects, therefore replacingR with The correction due to δR i is only relevant to NNLL order for one of the resolved emissions. This special emission is denoted by the subscript s in Eq. (A.26). After expanding to first order the corresponding term proportional to δR i arising from the initial R factor, one ends up with the second term in the curly bracket in Eq. (A.26), see Ref. [26] for a full derivation. The same manipulations apply to theR correction coming from the expansion (A.17) discussed above. Moreover, we introduced the following generalized luminosity coefficient and its NLL approximation We now report all the various ingredients entering the above formulas. The O(α s ) correction to the collinear coefficient functions reads ij are the LO Altarelli-Parisi splitting functions with β 0 = (11C A − 2n f )/(12π) and P (0), ij (z) are given by The functionH (1) where H (1) denotes the finite one-loop virtual correction to the gg → H process and A (i) , B (i) are reported below. For the top contribution in the m t → ∞ approximation, H (1) reads The result including full quark mass dependence has been computed analytically in Refs. [40,52,53]. 17 We expand the Sudakov radiator as We introduce and write (A.40) 17 In our implementation we take both the Born amplitude and the virtual corrections from Ref. [40].
The expressions ofR , δR , andR used in Eq. (A.26) are defined as The β function coefficients read Finally, we have

A.3 Matching to fixed order
In this section we discuss the matching of the resummed and the fixed-order results. We work at the level of the cumulative distribution Σ, that at NNLO reads (A. 45) We stress that in the main text we only show results for the differential p ⊥ distribution, therefore we label them as NLO. This corresponds to what we label as NNLO at the integrated level in this appendix. Since the total gg → H cross section is not known in the full SM beyond NLO, we approximate the NNLO correction to σ NNLO tot by multiplying the exact NLO result by the NNLO/NLO K factor as computed in the m t → ∞, m b → 0 limit. We stress, however, that at the level of the differential distributions we are interested in, this approximation is formally a N 3 LL effect, and it is beyond the accuracy considered in our study.
In order to assess the uncertainty associated with the matching procedure, we consider here two different matching schemes. The first scheme we introduce is the common additive scheme discussed in the main text defined as Σ add (p ⊥ ) = Σ NNLL (p ⊥ ) + Σ NNLO (p ⊥ ) − T NNLO Σ NNLL (p ⊥ ) . (A.46) Since the O(α 2 s ) (relative to the Born) collinear coefficient functions and virtual corrections are unknown in the full SM, in the additive scheme we approximate them by multiplying the HEFT ones by the exact Born squared amplitude.
The second scheme we consider belongs to the class of multiplicative schemes. In the text, we schematically defined it as Σ mult (p ⊥ ) = Σ NNLL (p ⊥ ) T NNLO Σ NNLO (p ⊥ ) Σ NNLL (p ⊥ ) .
(A. 47) We recall that we indicate with T NNLO [f ] the fixed-order expansion of f to NNLO. The two schemes (A.46), (A.47) are equivalent at the perturbative order we are working at, as they only differ by N 3 LO and N 3 LL terms. The main difference between the two schemes is that, in the multiplicative approach, unlike in the additive one, higher-order corrections are damped by the resummation factor Σ NNLL at low p ⊥ . One advantage of the multiplicative solution is that the NNLO constant terms, of formal accuracy N 3 LL, are automatically extracted from the fixed order in the procedure. Furthermore, as we explained in the text, in this case higher order effects introduced by the resummation follow the same scaling in p ⊥ of the fixed-order result, which at least partially mimics higher order form-factor effects. However, there is a drawback in using Eq. (A.47) as is. Indeed, Σ NNLL does not tend to one for p ⊥ Q, but rather to the luminosity factor defined in Eq. (A.28) evaluated atL = 0. Therefore, the fixed-order result Σ NNLO at large p ⊥ receives a relative spurious correction of order α 3 Despite being formally of higher order, these effects can be moderately sizable in processes with large K factors such as Higgs production. There are different possible solutions to this problem. In Ref. [26] the resummed factor (and the relative expansion) was modified by introducing a damping factor as Σ NNLL → Σ NNLL Z , (A. 49) where Z is a p ⊥ -dependent exponent that effectively acts as a smoothened Θ function that tends to zero at large p ⊥ . This solution, however, introduces new parameters that control the scaling of the damping factor Z (see Section 4.2 of Ref. [26] for details). In this article we adopt a simpler solution, which avoids the introduction of extra parameters in the matching scheme. We therefore define the multiplicative matching scheme by normalizing the resummed prefactor to its asymptotic value for atL → 0. This is simply given by This ensures that in the p ⊥ Q limit Eq. (A.51) reproduces by construction the fixed-order result, and no large spurious, higher-order, corrections arise in this region. The detailed matching formulas for the two schemes considered in our analysis are reported below.
We start by introducing a convenient notation for the perturbative expansion of the various ingredients. We define where, as above, we grouped the terms entering at NLO, and NNLO within curly brackets.