The Higgs transverse momentum distribution in gluon fusion as a multiscale problem

We consider Higgs production in gluon fusion and in particular the prediction of the Higgs transverse momentum distribution. We discuss the ambiguities affecting the matching procedure between fixed order matrix elements and the resummation to all orders of the terms enhanced by log(pTH/mH) factors. Following a recent proposal [1], we argue that the gluon fusion process, computed considering two active quark flavors, is a multiscale problem from the point of view of the resummation of the collinear singular terms. We perform an analysis at parton level of the collinear behavior of the Oαs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s\right) $$\end{document} real emission amplitudes; relying on the collinear singularities structure of the latter, we derive an upper limit to the range of transverse momenta where the collinear approximation is valid. This scale is then used as the value of the resummation scale in the analytic resummation framework or as the value of the h parameter in the POWHEG-BOX code. A variation of this scale can be used to generate an uncertainty band associated to the matching procedure. Finally, we provide a phenomenological analysis in the Standard Model, in the Two Higgs Doublet Model and in the Minimal Supersymmetric Standard Model. In the two latter cases, we provide an ansatz for the central value of the matching parameters not only for a Standard Model-like Higgs boson, but also for heavy scalars and in scenarios where the bottom quark may play the dominant role.


Introduction
A new state with a mass of approximately 125 GeV has been observed at the LHC [2,3]. Many investigations are under way to determine its properties and to test its compatibility with the Higgs scalar boson of the Standard Model (SM). The precise measurements of the total production cross section and of the branching ratios in the different allowed decay channels [4,5] have shown that the new state couples to the known fermions and gauge bosons following the SM predictions. Other studies target the kinematics of the decay products to distinguish among the various spin-parity combinations [6]. Finally, further work will be necessary to clarify the structure of the scalar potential.
In the SM the main production mode of the Higgs boson at hadron colliders is through the gluon fusion mechanism. The coupling of the Higgs boson to the gluons is mediated by a loop of colored particles, with the largest contribution to the process given by the JHEP01(2016)056 top quark. The gluon fusion cross section is very well approximated by a Heavy Quark Effective Field Theory (HQEFT), where the Higgs mass is considered very small with respect to the one of the top quark. The coupling of the Higgs boson to gluons is then proportional to the Fermi constant and to the strong coupling, but it is independent of the top Yukawa coupling. In this approach the very large NLO and NNLO-QCD corrections to the LO process (+100% and +30% of the LO result respectively) have been evaluated in refs. [7][8][9] and in refs. [10][11][12][13][14][15]. Recently, expressions for the N 3 LO corrections have been published in refs. [16][17][18][19][20][21]. The calculation using the complete SM Lagrangian was done up to NLO-QCD [22][23][24][25]. The exact treatment of the quark loops (mostly from top, bottom and charm) at NLO-QCD yields an O(−1%) correction, for an Higgs with a mass of m H 125 GeV and a collision energy √ S = 14 TeV. Moreover finite top-mass effects at NNLO-QCD have been estimated and found to be of O(1%) [26][27][28][29][30][31]. Beyond fixed-order QCD corrections, also soft-gluon resummation effects are available [32][33][34][35][36]. Moreover, the first-order electroweak (EW) contributions have been evaluated in refs. [37][38][39][40][41][42][43][44] and an estimate of the mixed QCD-EW contributions has been presented in ref. [45]. The PDF and α s uncertainties on the total Higgs production cross section have been studied in ref. [46].
The production cross section of a Higgs boson at large transverse momentum has been computed at LO-QCD, retaining the full quark-mass dependence, in refs. [47,48]. The NLO-EW corrections to this observable have been considered in refs. [49,50] in the HQEFT limit. The NLO-QCD corrections, in the HQEFT, have been computed [51][52][53]. An estimation of top-mass effects at NLO-QCD has been presented in ref. [54]. The first results towards the determination of the Higgs production at large transverse momentum, in the HQEFT, with NNLO-QCD accuracy, have been presented in ref. [55].
In this paper we want to reconsider the uncertainties that affect the theoretical prediction for the Higgs boson transverse momentum p H ⊥ . The transverse momentum distribution is an observable generated by QCD radiation. In the region of small p H ⊥ the presence of terms enhanced by large log(p H ⊥ /m H ) factors spoils the accuracy of the fixed-order results; in order to obtain a physically meaningful prediction these logarithms have to be resummed. Various techniques are available to perform the resummation. Once the latter is achieved, the resummed result has to be matched to the fixed-order one. Particular care is required to avoid the double counting of those logarithmic contributions that are present in both computations. The matching procedure introduces additional unphysical variables, the matching parameters, that define how the spectrum is divided into a soft region, where the resummed result is indeed applied, and a hard region where the fixed-order result is instead considered as the correct description of the spectrum.
In the HQEFT framework, the corrections up to NLO-QCD for the Higgs transverse momentum distribution have been analytically computed and matched with the transverse momentum resummation at NNLL accuracy. The results have been originally implemented in the code HqT [56][57][58] and later in the parton Monte Carlo program HRes [59]. A similar discussion, in the Soft Collinear Effective Theory (SCET) approach, has been presented in refs. [60][61][62]. In the context of matched NLO+Parton Shower (PS) Monte Carlo event generators, which implement the resummation algorithmically in the computer code, the results in the HQEFT, for Higgs production via gluon fusion, have been presented in JHEP01(2016)056 refs. [63,64]. Two shower Monte Carlo codes that retain the NNLO-QCD accuracy on the inclusive observables, in the HQEFT, have been presented in refs. [65,66].
Despite the fact that the exact matrix elements retaining the full dependence on the quark masses were available for quite some time, they have been implemented in a NLO+PS Monte Carlo for the first time in ref. [67], in the POWHEG approach, and later in MC@NLO [68]. A similar study, in the framework of analytic resummation, has been presented in ref. [69] and later in ref. [1]. Recently, these effects have been implemented in the NNLOPS code [70]. Quark mass effects have also been discussed for observables like the jet veto distribution in refs. [71,72]. Moreover, in ref. [71], the structure of the collinear singularities and of the regular terms present in gluon fusion at O(α s ) is analyzed in detail.
In ref. [1] it has been pointed out that the matched computation of the Higgs transverse momentum distribution is a problem with three scales, namely the Higgs mass, the internal quark mass and the transverse momentum of the Higgs boson. The matching prescription between fixed-order and resummed results should account for all these scales, to avoid, as far as possible, the inclusion of spuriously large higher-order terms in the final result. It should be noted that the presence of non-negligible interference effects between the top and the bottom quarks assigns a simultaneous active role to both internal quarks present in the scattering amplitude.
In the framework of SCET, the separation between the singular regions where a resummation is needed and the corresponding regular parts has been discussed in refs. [73,74] with the introduction of appropriate profile functions at the level of the hadronic cross section; this approach has been applied to Higgs studies in ref. [75]. The problem of the determination of a sensible value for the scale that separates the two transverse momentum regions, the one where the resummation is needed and the one where a fixed-order description is reliable, has been discussed in QCD, at the level of the partonic cross section, in ref. [76]. Recently, in ref. [77], the determination of these scales has been realized in QCD, with an approach that exploits some general properties of the Higgs transverse momentum distribution at hadron level, to derive the largest interval of transverse momenta where the resummed expression can be applied.
In this paper we elaborate the approach of ref. [76], and present a derivation at parton level of the interval of transverse momenta where the collinear approximation of the squared matrix element is accurate and the transverse momentum resummation can be safely applied. A comparison of the present results against those of refs. [77,78] is currently ongoing [79].
Higgs production via gluon fusion may provide interesting information about possible signals of physics beyond the Standard Model (BSM), like those possible in the Two Higgs Doublets Model (2HDM) or those predicted in the Minimal Supersymmetric Standard Model (MSSM), thanks to the possible exchange in the loop of new colored particle that act as mediators of the interaction between the gluons and the Higgs boson. The total cross section for Higgs production (see ref. [80] for a recent review) and the Higgs transverse momentum distributions provide complementary information (see e.g. ref. [67]) to disentangle the SM from MSSM. The possibility of extracting sensible information from the data depends on the accuracy of the prediction of the p H ⊥ distribution, and, among others, on the choice of the matching parameters.

JHEP01(2016)056
The outline of the paper is the following. In section 2 we recall the basic elements needed to formulate the transverse momentum resummation and to match the corresponding expression with fixed order results; we make some comments on the analytic procedure and discuss in more detail the NLO+PS Monte Carlo formulation. In particular we discuss in both cases the role of the scales associated to the matching and we describe the differences between the two approaches. In section 3 we discuss in detail the gg → gH and the qg → qH processes, with respect to their collinear behavior. The latter is used to identify an interval of transverse momenta where the collinear approximation of the squared matrix element is accurate and where it is thus safe to apply the resummation procedure; we introduce a scale w that represents the upper bound of this interval. We discuss both scalar and pseudoscalar Higgs boson production and determine, as a function of the Higgs and the quark masses, in a model independent way, the scale w, which constitutes our main result. In section 4 we perform a phenomenological analysis in the SM: the numerical results of the previous section are applied in the analytic resummation context, with the code HRes, and in the NLO+PS Monte Carlo framework, with the code gg H quark-mass-effects present in the POWHEG-BOX; the corresponding Higgs p H ⊥ distributions are eventually compared. Finally, in section 5 we discuss the implications of the determination of the scale w in the MSSM and in the 2HDM, with the possible production of new heavy states with masses of several hundred GeV and with a possible strong coupling of the Higgs to the bottom quark, enhanced with respect to the SM case. For this study we use the generators gg H MSSM and gg H 2HDM, also present in the POWHEG-BOX. 2 Remarks on the computation of the Higgs p H ⊥ distribution

Analytic resummation and the collinear limit
The Higgs boson acquires a transverse momentum p H ⊥ because of its recoil against QCD radiation. In fixed-order perturbation theory the emission of initial state massless partons yields, in the collinear limit, a logarithmic divergence of the Higgs transverse momentum distribution, signaling a breakdown of the perturbative approach, with an effective expansion parameter α s (p H ⊥ ) log(p H ⊥ /m H ) ∼ 1 in the phase space region of vanishing p H ⊥ . The analytic resummation to all orders of the terms α s (p H ⊥ ) log(p H ⊥ /m H ) n is performed by exploiting the universal properties of QCD radiation in the collinear limit and restores an acceptable physical behavior (the Sudakov suppression) of the Higgs transverse momentum distribution in the limit p H ⊥ → 0 [81][82][83][84][85][86][87][88][89][90]. In the collinear limit p H ⊥ → 0 the amplitude for the real emission processes gg → gH and qg → qH diverge and can be written, via a Laurent expansion, as M exact = M div /p H ⊥ + M reg . In this limit, the second term can be neglected with respect to the first one and it is possible to recognize that M div is proportional to the Born amplitude times the appropriate radiation term. This factorized structure of the amplitude, neglecting the contribution coming from M reg which is assumed to be small, can be extended to all orders and it forms the basis of the resummation procedure. Indeed, we can iterate this factorization in the case of the amplitude for the emission of n additional partons. In impact parameter space, this procedure leads to a factorized expression with n divergent emission factors times a JHEP01(2016)056 term proportional to the Born amplitude. The expression for the approximated amplitude describing the emission of up to n partons can be cast in the form of an exponential series, which can thus be summed to all orders. The relative contribution of M reg to the full amplitude can be used to assess the accuracy of the collinear approximation and of the factorization hypothesis.
The resummed partonic cross section has a factorized structure given by the product of a universal exponential factor, which accounts for the resummation to all orders of the logarithmically divergent terms, multiplied by a process dependent function, which describes the details of the hard scattering process. This factorization requires the introduction of a scale µ res , called resummation scale [56]. The latter defines the region where the resummation is applied and it is usually set to a value between 0 and the hard-scattering scale. A customary choice in the literature, for inclusive Higgs production, is to set the central valueμ res = m H /2 [56]. The precise choice of this value is one of the main topics of this paper and will be further discussed in the next sections. Analogously to what happens with the renormalization and factorization scales, the physical observables should not depend upon µ res , but the truncation at a fixed order of the logarithmic expansion leaves a residual dependence on it, which can be used to estimate the uncertainty due to the missing higher-order logarithmic terms; a variation of the scale µ res in the interval [μ res /2, 2μ res ] is customarily adopted.
The matching procedure requires to fix the integral of the Higgs transverse momentum distribution to a constant, which is conventionally set to the value of the fixed order total cross section [56]. This constraint holds exactly for any choice of µ res , so that any variation of the resummation scale modifies the shape of the distribution but not its integral and yields thus a correlation between low-and intermediate-p H ⊥ regions.

Numerical resummation in the NLO+PS framework
Another approach to the resummation of terms enhanced by the factor log(p H ⊥ /m H ) is the one obtained in the context of PS Monte Carlo, where the multiple emission of partons is numerically simulated via the PS algorithm. The matching between the fixed order NLO-QCD results and the PS has been discussed in refs. [63,91,92] and it is implemented in several tools regularly used in the experimental analyses.
In a sufficiently general way we can write the matching formula as The phase space is factorized into the product of the Born and the real emission components, dΦ = dΦ B dΦ r . The Born squared matrix element is denoted by B whileB is the NLO normalization factor. The latter is defined as In this formulaV fin represents the UV-and IR-regularized virtual contribution. We use the hat to indicate that an amplitude has been IR-regularized. The partonic subprocesses

JHEP01(2016)056
with the emission of an additional real parton can be split into two groups: those that are divergent in the limit of collinear emission, called R div , and the ones that are instead regular, R reg . We can further subdivide the squared matrix elements of the divergent subprocesses in two parts: The term R s contains the collinear singularity of R div , while R f is a finite remainder. Finally, we use the symbol ∆ s t for the Sudakov form factor, with t as the shower ordering variable: The splitting of R div in eq. (2.3) is defined up to a finite part which can be reabsorbed in R s . In the literature two different choices have been adopted: in POWHEG R s = R div , while in MC@NLO R s ∝ α s P ij B is proportional to the product of the Born matrix elements times the relevant Altarelli-Parisi splitting functions.
It is interesting to observe that different definitions for R s generate higher-order effects in the matched differential cross section. The possibility of defining the finite part R f in an arbitrary way can be exploited to parameterize the uncertainties related to the matching procedure.

The role of the damping factor D h in the POWHEG-BOX framework
In the POWHEG-BOX framework, the separation between R s and R f can be achieved in a dynamical way using the damping factor D h , defined as (2.5) The divergent and the regular part of R div = R s + R f are then defined as: The role of the scale h is to separate the low and the high transverse-momentum regions and it therefore specifies the range of momenta for which the Sudakov form factor is possibly different from 1. In the limit p H ⊥ h we obtain R s → R div and R f → 0. In this limit the Higgs p H ⊥ distribution is suppressed by the Sudakov form factor. On the other hand, when p H ⊥ h we have R s → 0 and R f → R div and the Sudakov form factor tends to 1. In this latter regime the emission of a real parton is described at fixed order by the matrix elements R f = R div .
The differential distribution generated according to eq. (2.1) contains higher order terms, beyond the claimed accuracy of the calculation, due to the product ofB × R s . Indeed in the large p H ⊥ region we have Originally the factor D h was introduced to damp the R s contribution at large p H ⊥ and to recover the exact fixed order result in this kinematic region, at the level of the first emission handled by POWHEG.
By varying the scale h, it is possible to check how well the fixed order distribution is recovered for large values of p H ⊥ , as can be seen from figure 1. We observe that, while at the level of the first emission generated by POWHEG (obtained at the level of Les Houches Event File (LHEF)) the NLO result is fully recovered, the showering of the events causes the high-p H ⊥ tail of the distribution to rise over the NLO prediction.
The total NLO cross section is always preserved for any value of h, as can be checked by integrating eq. (2.1) over the whole phase space. This property implies in turn that the low-and high-p H ⊥ regions of the differential cross section are correlated. Any increase of the distribution at low-p H ⊥ translates in a decrease of the high-p H ⊥ tail and vice versa.

JHEP01(2016)056
The role effectively played by the scale h has some similarities with the one described in section 2.1 for the resummation scale µ res : indeed, for p H ⊥ < h or for p H ⊥ < µ res the Sudakov suppression yields a regular behavior of the Higgs transverse momentum distribution, whereas for p H ⊥ larger than these scales the fixed-order description is recovered, at the level of description given by POWHEG. It should however be remarked that µ res and h have a completely different origin. The scale µ res is introduced as the scale at which the resummation is defined and the factorization of the partonic cross section implemented. It necessarily appears in the arguments of the logarithmic terms that are resummed. The damping factor D h is instead a convenient p H ⊥ -dependent parameterization of the ambiguity in the definition of R s . In a different perspective, the scale h controls the range of p H ⊥ over which the first term in eq. ( 2.1) is active in the generation of the first real emission. Since this term contains the normalization factorB, the scale h in turn controls also how the total NLO cross section is spread over the p H ⊥ distribution.

The value of the SCALUP variable
The emission of the radiation in the POWHEG approach is described by eq. (2.1). Neglecting the contribution coming from the term R reg (negligible in the case of the Higgs production in gluon fusion), we have two different categories of events. One corresponds to the terms in curly brackets (B-events), while the second one is described by the term R f (remnant events). The latter is present only if the damping factor D h is used.
To avoid double counting of the emissions, in the POWHEG approach the PS is required to emit radiation at transverse-momentum scale lower than the one of the parton emitted by the POWHEG-BOX. More in detail, in the case ofB events the PS should start to consider the possibility of an emission exactly at the scale at which the POWHEG parton was radiated. In the default POWHEG-BOX implementation the same choice is applied, for a reason of uniformity, also to the remnant events. This information is therefore computed on an eventby-event basis and the passed to the PS using the SCALUP field in the LHE event record.
It might happen that the value of SCALUP is large and that the description by the PS of real radiation at large transverse momenta is not accurate, since this approach is based on the soft/collinear approximation. It is then natural to consider as an option the possibility of setting an upper bound to the value of the SCALUP variable, close or equal to the one adopted for the h parameter. 1 This choice is applied only to the events generated by the R f part of the real matrix elements, in order to respect the POWHEG accuracy given by the first term in eq. (2.1) and account for a higher-order effect. Since these events are relevant only for the description of the high-p H ⊥ region (in turn defined by the scale choice h), we do not expect a modification of the shape of the distribution in the low-p H ⊥ region. Indeed, in this way the action of the PS is restricted to the lower part of the p H ⊥ spectrum, whereas the large-p H ⊥ tail is described purely by the LO matrix elements. In section 4 we will compare the description of the Higgs high-p H ⊥ tail with the default and with the modified SCALUP values, in the case of the SM.

JHEP01(2016)056 3 Collinear approximation of partonic squared matrix elements
In the previous section we have recalled that the resummation to all orders of the terms enhanced by a log(p H ⊥ /m H ) factor is possible thanks to the factorization of the squared matrix element in the collinear limit p H ⊥ → 0. Based upon those considerations, we now explain our procedure to determine the accuracy of the collinear approximation of the full squared matrix element with respect to p H ⊥ , focusing for simplicity on the channel gg → gH. We then derive numerically the value w of the upper limit of the p H ⊥ range where the collinear approximation is accurate, for the scalar and the pseudoscalar final states, considering both the gg → gH and the qg → qH channels.

Helicity amplitudes and kinematic variables
We consider the helicity amplitudes 2 M λ 1 ,λ 2 ,λ 3 (s, p H ⊥ , m 2 H ) for the process gg → gH, whose complete expressions can be found for example in ref. [48]. We reorganize them, via a Laurent expansion, as follows: and we use this decomposition to compute the unpolarized squared matrix element exactly, |M| 2 , or its collinearly divergent part |M div /p H ⊥ | 2 . We define the ratio C: which quantifies how the unpolarized exact squared matrix element differs from its collinear approximation as a function of p H ⊥ . We observe that by construction we have lim p H ⊥ →0 C(s, p H ⊥ , m 2 H ) = 1. In our study we also consider the behavior of the interference term between the top and the bottom quark. For this specific case we redefine the parameter C as We introduce the following practical criterion: the regular part of the amplitude becomes non-negligible with respect to its collinear counterpart for a value w of p H ⊥ such that To fix the setup of our study we chooseC = 0.1. This value is arbitrary, but its order of magnitude can be justified in the framework of a QCD calculation, since the size of the terms without a collinear logarithmic enhancement is α s /π times a coefficient of order 1. We do not assign any special meaning to the scale that will be found with our analysis but we rather consider it as a starting point to compute an uncertainty band. At the end of the section we analyze the dependence on the specific value ofC.

JHEP01(2016)056
The amplitude of the process gg → Hg is a function of two independent kinematic variables, e.g. s and p H ⊥ . The production of a final state with a definite p H ⊥ requires a minimum value for s: We study the behavior of the amplitude as a function of p H ⊥ for s = s min + s soft , where s soft is necessary to avoid the soft divergence and focus only on the collinear behavior. The choice of a value of s close to s min is phenomenologically motivated by the strong PDF suppression in the hadronic cross section for increasing partonic s.
An analogous procedure is used to determine the scale w for the qg → qH subprocess, with the analytic expressions of ref. [25] and in the case of pseudoscalar production using the formulae in ref. [93].

Partonic analysis
We assume that the full amplitude is the sum of a top and a bottom contribution, neglecting the light quark generations. Furthermore, we do not consider the possibility of additional colored particles running in the loop, since the current LHC results hint to the fact that, if these states exist, their mass is probably much larger than the top mass. Therefore they would not affect the shape of the Higgs transverse momentum distribution for values of p H ⊥ that are phenomenologically interesting.
Under the assumption that the coupling of the gluons to the quarks is the one dictated by QCD and that all the details about the coupling of the Higgs to the quarks can be factorized from the rest of the amplitude, we can consider the value of the scales w t , w b and w i , that we respectively find in the case of squared matrix elements with only top quark diagrams, only bottom diagrams or for the top-bottom interference, as model independent. As a consequence, the determination of the scales depends only on the quark and the Higgs masses.
While the scales computed with only one quark might have a physical interpretation in the BSM scenarios where that quark yields the dominant contribution to the cross section, the scale of the interference terms, while unphysical, is a necessary tool to treat accurately the full theory in scenarios where both top and bottom quarks are equally important. Since the full squared matrix element, including top and bottom quarks, factorizes in the collinear limit, the same pattern should be followed not only by the terms with the squared amplitude of one single quark, but also by the interference terms, making our treatment viable.
In sections 4 and 5 we will discuss how these results can be exploited in a model specific framework.

Scalar Higgs
To exemplify the outcomes of our procedure, we show the results for the variable C We first discuss the impact of the regular terms in the case of a light Higgs. We compare the results obtained with the exact matrix elements including only the top quark with the ones in the HQEFT; we observe that in both models a deviation by more than 10% from the collinear approximation occurs for p H ⊥ > 55 GeV. Since it is present in both cases, this effect should thus not be interpreted as a top mass effect; the latter becomes visible for p H ⊥ > 150 GeV. From the analysis of the helicity amplitudes, we observe that this deviation from the collinear approximation stems from M −+− . For the bottom quark, the deviation from the collinear approximation starts from p H ⊥ > 19 GeV. In the case of the interference terms, we observe that the determination of the scale w i is dominated by the behavior of the bottom amplitude; the corresponding value, w i = 9 GeV, is smaller than the ones obtained in the other two cases.
In the case of a heavy Higgs, with m H > m t , m b , the scale of the process is set by the mass of the boson (e.g. m H = 500 GeV) and the HQEFT approximation of the amplitude is not valid. We observe that the amplitude that includes only the top-quark diagrams deviates from its collinear approximation 3 for p H ⊥ > 111 GeV. Instead, the squared matrix element that includes only the bottom-quark diagrams deviates from its collinear approximation for p H ⊥ > 63 GeV. Finally, for the interference terms we find the bound p H ⊥ > 18 GeV.  125  55  19  9  24  7  5  48  18  9  200  85  29  16  21  5  5  71  27  14  300  132 41  25  17  4  4  111  38  23  350  102 47  28  15  4  4  87  43  26  400  94  52  26  14  4  3  81  49  23  500  111 63  18  13  3  2  96  58  17  600  133 73  6  13  3  0  113  68  6   the bottom diagrams or for the interference of the top and bottom amplitudes; the results are presented separately for the two partonic subprocesses, gg → gH and qg → qH.
We observe that in the gg → gH channel both scales w t,b increase with the Higgs mass, with the exception of the region of real top-pair production threshold, where the effect on w t of additional terms that induce a deviation from the collinear approximation is visible. In the bottom-quark case such phenomenon does not show up, because for realistic values of m H the process scale is always well above the bottom-pair production threshold. The interference scale w i has a peculiar behavior: in fact, it shows a growth with m H until the top-pair production threshold and then it decreases for larger m H , until it vanishes for m H = 589 GeV, with our m t and m b choices; for even larger m H values it grows JHEP01(2016)056 again. In order to explain why w i vanishes, we should recall that the interference terms, as a function of m H and for fixed m t and m b , are not positive definite and may change sign for a specific value of m H ; in particular, when the underlying LO (i.e. of the process gg → H) interference terms vanish, also the collinear approximation does. In this point the interference terms of the processes gg → gH and qg → qH are thus collinear finite, the function C(s, p H ⊥ , m H ) diverges for all p H ⊥ and the scale w i is equal to zero, indicating that the p H ⊥ distribution is regular and a LL resummation is not needed. It should be noted that, for this specific configuration, the importance of the interference term is in any case small, since it vanishes at LO.
We observe that in the qg → qH channel the scales are lower than in the previous case and that they decrease for increasing values of the Higgs mass. The basic argument to explain this different behavior can be found analytically in the HQEFT: we expand the ratio C(s, p H ⊥ , m 2 H ) in powers of p H ⊥ around p H ⊥ = 0 and we find The different behavior with respect to m H of the scale w is due, in the gg case, to the fact that the function C receives corrections with negative powers of m H , so that for heavy Higgs masses there is a larger interval of p H ⊥ where the collinear limit provides a good approximation of the full result; in the qg case instead, there are corrections quadratic in m H , such that the deviation of C from 1, for large m H , occurs at smaller p H ⊥ values. A numerical analysis with the full dependence on the top and bottom masses confirms the explanation derived above in the HQEFT.
The interference scale vanishes, as expected, for the same m H value in the gg → gH and the qg → qH channel, since they factorize to the same LO term.
The different values of the scales w t,b obtained in the two partonic channels gg and qg give rise to a practical problem, in case one wants to use at hadron level one single scale to control the effects of multiple parton emissions; given that the w value from the gg channel is always larger than the one from the qg channel, we can expect that the final value will lie in between; we evaluate it with a weighted average, with the relative contributions of the two channels in each bin, further adjusted to account for the shape of the physical distribution. We define  125  60  19  11  24  7  6  52  18  10  200  126 29  18  22  5  5  102  27  16  300  122 41  28  18  4  4  103  38  25  350  82  47  25  15  4  4  70  43  23  400  99  52  15  14  4  2  86  49  14  500  127 63  15  12  3  2  109  58  14  600  155 73  36  11  3  51  132  68  In figure 3, and in table 1, in the last three columns to the right, we show the results of this combination, which are our best determination for the scales to be used in the simulation of the hadronic differential cross section. 4 We have used the code SusHi [94], with √ S = 13 TeV, to compute the weights used in eq. (3.8). Since this procedure requires the evaluation of the hadronic cross section, the combined scales are dependent on the √ S value used and on the other hadronic parameters. In particular this is true for the choice of the renormalization and factorization scale, that we have assumed to be µ r = µ f = m H . However we have verified that the effect on the channel-combined value for the scales is only at the of few GeVs, well within the uncertainty band that we are considering. A finer scan in the Higgs mass, is provided in tables 4 and 5, in appendix A.

Pseudoscalar Higgs
In table 2 we present a sample of the results, analogous to the ones of the previous subsection, for the case of pseudoscalar Higgs production.
The general behavior of the two partonic channels is similar to the one observed for scalar production. One difference can be observed at the top-pair threshold, where a cusp appears in the w t prediction, reflecting the analogous feature of the total cross section. The scale w i vanishes for a different value of the pseudoscalar mass, m A = 445 GeV, because of the different LO dependence on m A , m t and m b . As for the scalar case, a more detailed scan as a function of m A is available in tables 4 and 5, in appendix A.

Dependence on auxiliary parameters
The value of the resummation scale has been determined with an analysis of the partonic squared matrix element, for fixed value of the partonic invariant s. For a given final state configuration and in particular for a given value of p H ⊥ , the hadronic distribution receives contributions from all the partonic cross sections with s min ≤ s ≤ S, where S is the 4 The relative weight of the two partonic channels, as a function of the Higgs mass, is slowly varying, so that we can approximate the result of equation ( Due to the assumptions used in our procedure, we stress that the determination of the central value for w does not have an absolute meaning. It is rather the starting point to define an interval of reasonable values for the scale w that in turn should be used to compute an uncertainty band for the transverse momentum distribution.

Standard Model phenomenology
We consider now the evaluation of the Higgs transverse momentum distribution in protonproton collisions at the LHC in the SM. We use the analytic results of [95] implemented in the public code HRes and the shower Monte Carlo implemented in the POWHEG-BOX [67]. For the former, we study the impact of different choices of the resummation scale µ res ,

JHEP01(2016)056
while with the latter we vary the value h which enters the damping factor D h . In both cases we consider the possibility of a separate treatment of the top and of the bottom quark contributions. In the numerical analysis we use m t = 172.5 GeV, m b = 4.75 GeV, the PDF sets MSTW2008nlo68cl and MSTW2008nnlo68cl [96] with their corresponding values of α s (m Z ). We chose µ R = µ F = m H as the renormalization and factorization scales. We use PYTHIA8 [97,98] with the tune AU-CT10 to shower the POWHEG events. This specific tune was chosen since it is the same used by the ATLAS collaboration for their Higgs analyses. The center of mass energy at the LHC has been assumed to be √ S = 13 TeV.

Comparison of POWHEG and HRes
The gluon fusion process, including the top and the bottom quark diagrams in the scattering amplitude, is a three-scale problem, as was already stressed in ref. [1] and as we have seen in the previous sections: the Higgs mass, the value of p H ⊥ and the mass of the quark. The bottom quark contributions spoil the validity of the factorization hypothesis for p H ⊥ values smaller than in the top quark case and require a dedicated treatment. In order to make explicit the role of the top and of the bottom quarks, the squared matrix elements can be rearranged as where we have put in round bracket the quarks that run in the loops of the diagrams. The square brackets contain the top-bottom interference terms and the square of the modulus of the bottom amplitude. The rationale behind this rearrangement is that in the SM the dominant contribution to the gluon fusion is due to the top quark diagrams, while the bottom quark diagrams yield a correction to the former; it is thus reasonable to make one dedicated scale choice for the top quark and a second scale choice for all the other terms, even if they still include top quark diagrams via interference terms. We recall that by construction the total cross section does not depend on the value of the resummation scale in HRes (or equivalently of the scale h in POWHEG). This fact allows us to write the following identity where here and after, with a slight abuse of notation, we have introduced the symbol σ(q, µ) to indicate the total cross section evaluated with the quark q in the loops, using, in the numerical code, the matching parameter at the scale µ. The latter is the resummation scale Q i in HRes and the scale h in POWHEG. This equation is trivial for the total cross section, and represents a possible recipe for the evaluation of differential observables, specifically the Higgs boson transverse momentum. 5 For our phenomenological analysis we use two scales, one for the squared matrix element with only the top quark and one for the other contributions, to allow a comparison 5 In POWHEG, at the differential level, the extraction of a specific contribution by subtraction is bound to introduce spurious terms due to the fact that the Sudakov form factor is non-universal. However, due to our specific scale choices that guarantee a good accuracy of the collinear approximation in the p H ⊥ range where the Sudakov form factor has its major effect, we can argue that in this region the argument in the exponent of the Sudakov factor is well approximated by the relevant universal expression R/B αsPij/t, limiting the impact of the spurious terms.

JHEP01(2016)056
with the results presented in ref. [1]; we use a combination analogous to the one of eq. (4.2) to evaluate also the differential distributions.
In section 3.3 we have given an estimation of the uncertainty in the determination of the scales w by varying the auxiliary parameters that we have used in our computation. In theory it is possible to use the range of scales obtained with such a procedure as the range of values to be used for the matching parameter to estimate the uncertainty on the prediction for the transverse momentum distribution. However these values depend in a non-trivial manner on the Higgs mass considered. We observe that a variation by a factor of 2 of the central value widely covers the range of scales that we find with our explicit computation, thus yielding a conservative assessment of the uncertainty. To simplify the uncertaintyestimation procedure we have then decided to compute the uncertainty bands using the following standard prescription: we consider the 9 combinations of the pairs (µ t , µ b ) of the two matching parameters, which can be obtained from the sets (μ t /2,μ t , 2μ t ) and (μ b /2,μ b , 2μ b ), where we calledμ t andμ b the respective central values, and we take the envelope of all the predictions.
We consider the three following cases and, in each of them, we compute the uncertainty band according to the rule described above: The distributions obtained with HRes and POWHEG share the same matrix elements that describe at NLO-QCD the inclusive Higgs boson production, and differ by subleading NNLO and by higher order terms, which might nevertheless be numerically relevant.
The comparison of the shape 6 of the p H ⊥ distribution, in figure 5, of the results of item 1 (blue dot line) and 3 (dashed red line) is meant to expose the differences of the two codes taken with their default setup, when they are run with the same accuracy for the total cross section, NLO-QCD, and with the same value for the matching parameters. On the left we show the absolute comparison of the results, while on the right we show the ratio of the different predictions over the one obtained with POWHEG and the HRes scale choice (item 1).
As discussed in section 2 the two basic formulae used to generate the Higgs p H ⊥ spectrum differ by subleading O(α 2 s ) and higher-order terms, part of which are controlled by the resummation scale in HRes or by the h scale in POWHEG. For the above reason, even if we assign the same numerical values to the scales Q and h, we expect a certain level of discrepancy for the central predictions. Indeed we see that in the region where resummation effects are relevant, the two codes behave differently, with HRes giving a softer distribution than POWHEG. Specifically, the shape of the distribution produced by HRes is larger than the one from POWHEG for p H ⊥ ≤ 50 GeV, while for higher p H ⊥ the behavior is the opposite. In the high-p H ⊥ region, for p H ⊥ ≥ m H , we see that the HRes result coincides with the fixed-order distribution: in fact, the code HRes uses the full matched expression for p H ⊥ values smaller than m H and implements a smooth transition to the pure fixed-order expression, which is used in the high-p H ⊥ tail; for this same reason, the HRes resummation scale uncertainty band vanishes in this part of the spectrum.
In the high-p H ⊥ range POWHEG shows a distribution harder than the fixed-order one, because of the showering effects applied on top of the POWHEG formula for the first emission.
Since HRes does not include non-perturbative effects, which are present in the selected tune of the PYTHIA shower, an additional problem in the comparison emerges: the nonperturbative effects are relevant at small transverse momenta of the radiated partons. In addition, in the low-p H ⊥ region, the different expression of the HRes and POWHEG Sudakov form factors (for the latter see equation (2.4)) has a role to determine the precise shape of the distribution. By construction, the unitarity constraint, that forces the total cross section to be always preserved, implies an anti-correlation between the low-p H ⊥ and the high-p H ⊥ parts of the spectrum. The comparison in figure 5  In figure 6 we present the impact in POWHEG of a different choice of the variable SCALUP, as discussed in section 2.3. We set SCALUP=h t , a constant value, while we keep unchanged all the other parameters and in particular the value of the scales h t,b in the damping factor D(h). The choice for the SCALUP value is in accordance with the dominant role played by the top-quark loop in the SM. We observe that the central prediction of this modified POWHEG version is lower than the default one for p H ⊥ ≥ 200GeV and tends to recover the fixed-order distribution at large transverse momenta. We interpret the reduction of the differential cross section at large p H ⊥ as due to the missing contribution in this region from the PS emissions. The accuracy of the latter is questionable, since the PS is based on the soft/collinear approximation and might be inadequate to describe large-p H ⊥ radiation.

Beyond SM phenomenology
The description of the Higgs transverse momentum distribution in the SM, with m H = 125 GeV, is characterized by the dominant role played by the top-quark contribution, such that the bottom-quark effects can be treated as a correction. Moreover, with a light scalar Higgs, the HQEFT limit is a good approximation of the full SM, and the determination of

JHEP01(2016)056
the scale of validity of the collinear approximation (and hence of the applicability range of the resummation techniques) reduces to a problem involving only m H and p H ⊥ . At variance with the previous case, and still in the SM, we know that with a heavy Higgs boson, the description of the p H ⊥ distribution is a multiscale problem; indeed, the minimal energy scale necessary to produce the final state immediately probes the top-quark loop.
In a generic BSM scenario it is possible to consider enhanced couplings of the bottom quarks to a relatively heavy Higgs boson, scalar or pseudoscalar. In these configurations, our intuition, accustomed to a light SM-like Higgs phenomenology, may fail in the determination of the correct regime where the resummation techniques can be safely applied. Since a priori we do not know exactly how the contributions from the different quarks interplay in the full result, following ref. [77], we can generalize eq. (4.2) to where the last term allows us to use a separate scale for the top-bottom interference term. As before, the parton level analyses discussed in section 3 provide a model independent ansatz for the three relevant scales, µ t,b,i : these are the scales w t ,w b and w i , listed in tables 1 and 2 as a function of the Higgs boson mass. In order to illustrate the phenomenological consequences of our study, we show our predictions in the 2HDM and in the MSSM 7 and we compute the uncertainty bands with an extension of the procedure described in section 4: we consider all the 27 combinations of the three matching scales and then take their envelope. The range of scales spanned represents again a conservative choice to assess the matching uncertainty.

2HDM phenomenology
We consider the type-II 2HDM. We adopt a purely heuristic approach to show the impact of our study, choosing the 2HDM parameters that are relevant for the gluon fusion process by following only the requirement that they represent three different scenarios: one where the cross section is dominated by the top-quark; one where the contribution of the top and the bottom quark are of the same order of magnitude; and one where the process is dominated by the bottom quark matrix elements. The explicit values for the parameters are reported in table 3 for all the three scenarios. In all three cases we choose to study a heavy Higgs of m H = 500 GeV. The corresponding values for the scales are w t = 96 GeV, w b = 58 GeV and w i = 17 GeV. For the simulation we adopt the Monte Carlo generator gg H 2HDM available in the POWHEG-BOX.
We now present our best predictions obtained with the three-scale procedure and check how well they are approximated by a one-scale approach.
In figure 7 we show the results for the first scenario. In this case we have that the process is dominated by top quark contribution. Indeed we notice that the three scales  On the other hand, in the second case shown in figure 8, we have that the contributions coming from the two quarks are of the same order of magnitude. In this case we observe that the result obtained by using three scales is not recovered by simulations with just a single scale, with either the top or the bottom one.
Finally, in figure 9 we see that in the bottom dominated scenario, we have a similar situation as in the top dominated case, though the scale to be used in a one scale run is w b instead of w t .
In all three cases we stress that the using values of the order of m H /2 or m H /1.2 for the matching parameter h yields results that are in the best case at the limit of the uncertainty band.

MSSM phenomenology
We consider an explicit example in the MSSM by taking, in its parameter space, a point still allowed by the most recent available data, according to the analysis of ref. [99], and to the results of the code HiggsBounds [100][101][102][103]. The same point has been considered also in ref. [77]. We choose the so called m mod+ h scenario defined in ref. [99] and set M A = 500 GeV and tan β = 17 to fully specify our input parameters; as a result we obtain that the masses of the two CP-even Higgses are respectively m h = 125.6 GeV and m H = 499.9 GeV. The corresponding values of the w scales are: w t = 96 (109) GeV, w b = 58 (58) GeV and JHEP01(2016)056 w i = 17 (14) GeV, for a scalar (pseudoscalar) boson. We use these values to set the µ t,b,i parameters that enter eq. (5.1). For the simulation we adopt the Monte Carlo generator gg H MSSM available in the POWHEG-BOX. In the simulation we include the full particle content of the MSSM. We do not expect an important contribution from the squarks because in this point of the MSSM parameter space their masses are, respectively, mt 1 = 876 GeV, mt 2 = 1134 GeV, mb 1 = 1007 GeV and mb 2 = 999 GeV. With this specific parameter choice, the light CP-even Higgs is similar to the SM scalar, not only for the total cross section, but also for the shape of the p H ⊥ distribution. The heavy CP-even scalar and the pseudoscalar bosons have instead different properties, because of the different coupling strength to the top and to the bottom quarks. In

Conclusions
The study of the Higgs transverse momentum distribution may provide important insights about the properties of the recently discovered scalar resonance. The theoretical prediction of this observable requires, in the region of small p H ⊥ values, the resummation to all orders of terms enhanced by powers of log(p H ⊥ /m H ), while at large values of p H ⊥ , fixed-order calculations provide the most accurate description available. The consistent matching of the two approaches requires the introduction of a momentum scale, that separates the soft and the hard p H ⊥ regions. Since the validity of the resummation formalism relies on the collinear factorization of the squared matrix elements describing real parton emissions, we investigated the accuracy of the collinear approximation in the gluon fusion process, in the presence of an exact description of the top and bottom quarks running in the virtual loop. The discussion involves three scales, namely the Higgs mass, the Higgs transverse momentum and the quark masses.
Relying on the collinear singularities structure of the O(α s ) real matrix elements, we determined, in a model independent way, as a function only of the Higgs and the quark parton-level analysis and can be eventually used in any hadron-level computation (analytic or Monte Carlo) of the Higgs p H ⊥ distribution, following eq. (5.1). They indicate the upper limit of the p H ⊥ range where the resummed part of the cross section can be evaluated in a reliable way, because of the good accuracy of the collinear approximation of the full squared matrix elements. They represent an ansatz for the matching scales, whose values do not have an absolute meaning, but are rather the starting points to build an uncertainty band.
The procedure to compute an uncertainty band is described in section 4 and offers a simple but quite conservative recipe to derive this band. A more aggressive approach would exploit the scales obtained with a variation of the parameterC ∈ [0.05, 0.2], as discussed in section 3.
Our analysis is relevant for an accurate prediction of the Higgs p H ⊥ distribution, both in the SM and in BSM scenarios. In the latter case, our approach allows us to decompose the different contributions to the p H ⊥ distribution, also in the presence of a non trivial interplay between the Higgs transverse momentum and the Higgs, top and bottom masses, for any generic ratio between the strength of the couplings of the Higgs boson to the top and to the bottom quarks.
The description of the Higgs transverse momentum distribution, based on the use of three different scales for the matching parameter, represents our best ansatz for this observable. We remark, however, that in various cases this result can be accurately approximated with only one run that uses one single scale, the one associated to the dominant contribution to the scattering amplitude. This conclusion is obviously possible only a posteriori.
We stress the impact of the matching scale determination with one final comment, relevant in the context of the searches for new heavy scalars, referring to the results shown in figures 10 and 11. Our procedure defines the scales w t,b,i , whose variation in a given range is then exploited to compute an uncertainty band of the distribution. The results presented in section 5 are obtained with a conservative choice for the range of scale variation, as described at the beginning of the same section. The use of a single-scale simulation, with the matching scale set equal to the commonly adopted SM value m H /2, can lead to predictions that lie outside of this most conservative uncertainty band described above and that can differ with respect to our best central value by 30-40% both in the low-and in the high-p H ⊥ tails of the distribution.

JHEP01(2016)056
A Scan over the Higgs mass of the scales w t,b,i In the appendix we include two tables with the values of the combined gg-qg collineardeviation scales, for scalar and pseudoscalar masses between 100 GeV and 500 GeV, separately for the top, the bottom and the interference contribution. The top pole mass has been set to 172.5 GeV, while the bottom pole mass is equal 4.75 GeV, following the prescription by the Higgs Cross section Working Group (HXSWG). The merging of the scales was implemented by using the information on the relative importance of the two partonic subprocess as given by the code SusHi [94]. The latter was run at √ S = 13 TeV using the MSTW2008nlo68cl PDF set and setting µ r = µ f = m H .  Table 4. Values of the scales w gg+qg t,b,i as a function of the scalar and pseudoscalar Higgs mass.

JHEP01(2016)056
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.