Transverse momentum resummation for Higgs production via gluon fusion in the MSSM

The resummed transverse momentum distribution of supersymmetric Higgs bosons produced through gluon fusion at NLO+NLL is presented, including the exact quark and squark mass dependences. Considering various MSSM scenarios, we compare our results to previous ones within the POWHEG approach. We analyze the impact of the bottom loop which becomes the dominant contribution to the gluon fusion cross section for a wide range of the parameter space for the pseudo-scalar and heavy Higgs.


Introduction
After the observation of a Higgs boson of mass 125 GeV [1,2], the measurement of its properties has become one of the central targets of the LHC. From the theoretical side, precise predictions for the production and decay rates of such a particle in various models are crucial to pin down its nature. An enormous effort has already gone into precision calculations of the total cross section as well as kinematical distributions (see Refs. [3][4][5] for an overview).
In the Standard Model (SM), the Higgs is predominantly produced via gluon fusion, where the Higgs-gluon coupling is mediated by a quark loop. Its cross section is about an order of magnitude higher than the sum of all other processes, which retain their importance through their additional final state particles and/or their specific kinematics. The gluon fusion process has been studied in great detail over the recent years, leading to a significant decrease of the related theoretical uncertainties. In particular, the use of an effective theory approach for the calculation of higher order corrections allows, loosely speaking, to determine the cross section one perturbative order higher than in the full theory. Within this approach, also known as the heavy-top limit, the top quark is assumed to be infinitely heavy. The total cross section in this approximation is known up to next-to-next-to-leading order (NNLO) [6][7][8] and even parts of the next-to-NNLO are known [9][10][11][12]. Electro-weak corrections and further effects beyond NNLO have been evaluated in Refs. [13][14][15][16][17][18][19][20][21] for example. The uncertainty induced by the heavy-top limit has been shown to be below 1% for the total rate at NNLO [22][23][24][25][26].
While the effects of the four lightest quarks as mediators of the gluon-Higgs coupling is negligible ( 1%) and therefore usually omitted, the bottom quark contributes at the 5 − 10% level to the total cross section at next-to-leading order (NLO) [27,28]. Due to the smallness of the bottom-quark mass, one cannot apply the same approximation as for the top-quark contributions to evaluate radiative corrections for the bottom loop, but typically includes the full quark mass dependence in the calculation.
Kinematical distributions of the Higgs boson provide an important handle on the determination of Higgs properties (see, e.g., Refs. [29][30][31]). One of the most important differential observables in this respect is the transverse momentum (p T ) distribution of the Higgs. Once sufficient statistics have been collected at the LHC, the comparison of the experimental result for this spectrum to its theoretical prediction in various models should allow for further restrictions of the allowed parameter space of these models.
In the SM, the NNLO transverse momentum distribution of the Higgs produced via gluon fusion at p T > 0 has already been known for some time in the heavy-top limit [32][33][34]. 1 Sub-leading top-mass effects have been considered in Refs. [35,36]. Furthermore, also the fully-differential cross section has been determined up to NNLO [37][38][39]. However, it is well known that those perturbative calculations break down for small transverse momenta due to the occurrence of logarithmically enhanced terms in p T . Only a resummation of these terms to all orders in α s provides a proper theoretical prediction at small values of p T .
Transverse momentum resummation at leading logarithmic (LL) and next-to-LL (NLL) accuracy for the gluon fusion process in the heavy-top approximation has already been performed a long time ago [40][41][42]. Ref. [43] introduced a matching procedure to consistently combine the resummed distribution and the fixed order cross section valid at large p T . Its application to the p T spectrum of the Higgs at NNLO+NNLL was implemented using the effective theory approach into the publicly available program HqT [43][44][45]. 2 Finite top-and bottom-mass effects on the resummed p T spectrum have been considered in the POWHEG [48] approach [49] and by analytic resummation through NLO+NLL [47,50,51]. The small bottom-quark mass m b introduces an additional uncertainty because terms ln(m b /p T ) appear at the amplitude level 3 , which are potentially large and could spoil the collinear and soft approximation already at p T m b [47]. The small bottom Yukawa coupling in the SM prevents this uncertainty from becoming too severe though.
Supersymmetric extensions are among the most popular theories beyond the SM. The minimal supersymmetric SM (MSSM) contains two Higgs doublets, which lead to five physical Higgs bosons, three of which are neutral and two charged. The production of the neutral MSSM Higgs bosons is typically dominated by either of two processes, gluon fusion or bottom-quark annihilation. For the theoretical status of the latter process, we refer the reader to Ref. [52], where the resummed p T distribution through NNLO+NNLL was obtained, and the references therein. 4 In this paper, we focus on the gluon fusion process in the MSSM, but our calculation will be applicable also in a general 2-Higgs doublet model (2HDM). The total Higgs production cross section in gluon fusion in the MSSM has been calculated up to NLO within the MSSM [49,[57][58][59][60][61][62][63][64][65][66][67]. The currently most accurate total cross section in the MSSM can be obtained with the publicly available program SusHi [68]. It is obvious that our approach is also directly applicable to all neutral Higgs bosons of the 2HDM. 5 Our goal is the analytically resummed p T spectrum of all three neutral MSSM Higgs bosons produced through gluon fusion at NLO+NLL. Since the bottom Yukawa coupling can be significantly enhanced with respect to the SM, the issue of a proper treatment of 2 A Monte Carlo approach, based on the same resummation formalism, to add resummation effects to the differential NNLO cross section with respect to the Higgs and its decay products was implemented into the program HRes [46,47]. 3 However, in the limit pT → 0, these terms ∼ ln(m b /pT ) vanish, and therefore collinear and soft factorization is preserved. 4 The NNLO pT distribution for bottom-quark annihilation is already known for a while [53][54][55][56]. 5 Concerning the total cross section in 2HDMs, see Ref. [69].
bottom-quark induced effects on the cross section becomes more important. We propose a pragmatic way to separately set the resummation scale of these terms and to derive an estimate of the residual uncertainty.
We compare our results to the ones of a similar earlier study [49], which calculated the transverse momentum spectrum within the POWHEG approach [70] in combination with a parton shower.
The paper is organized as follows: In Section 2, we briefly review the formalism for the resummation of contributions at small transverse momenta in the gluon fusion process and discuss the required theoretical quantities. Our procedure for choosing the resummation scale is described in Section 3. Section 4 lists the input parameters and defines a set of MSSM parameter points which we use for our analysis. It also describes the way we determine the theoretical uncertainties. Numerical results are presented in Section 5, where we analyze the p T spectra for all three neutral MSSM Higgs bosons in specific scenarios and the impact of the relative contributions ordered by the respective Yukawa couplings that enter the cross section. Section 6 contains our conclusions.

Resummation and Matching
Consider the transverse momentum (p T ) distribution of a color-neutral heavy particle of mass M produced via a 2 → 1 process in QCD. For p T M , a fixed-order expansion of the cross section in the strong coupling α s can be applied. In the limit p T → 0, however, large logarithms ln(p T /M ) appear at fixed order, which spoil the validity of the perturbative expansion. A proper prediction of the distribution at p T M can be obtained by resumming these logarithms to all orders in α s . Following Ref. [43], we split the p T -dependent cross section as dσ dp 2 where the resummed logarithmic contributions in p T are contained in the first term on the r.h.s., while the second term remains finite as p T → 0. Working at finite orders, the cross section can be cast into the following form: where "f.o."(=fixed order) denotes the perturbative, and "l.a." the logarithmic accuracy (to be defined below) under consideration. The imposed matching condition defines the logarithmic accuracy needed at a specific perturbative order, and vice versa. In Eq. (2), all terms ∼ δ(p T ) are contained in the first term of the r.h.s.; in practical calculations, one can therefore disregard such terms in the second and third term since they will cancel among each other.
The matching procedure as proposed in Ref. [43] induces a unitarity constraint on the matched cross section which implies that the integral over p 2 T reproduces the total cross section σ tot at fixed order: In the next section, we will address the evaluation of dσ (res) /dp 2 T .

The resummed cross section
The resummation of large logarithmic contributions is performed in the impact parameter or b space, given by the Fourier transform w.r.t. the transverse momentum: 6 [72,73] dσ F,(res) dp 2 with q ∈ {u, d, s, c, b}, a numerical constant 7 b 0 = 2 exp(−γ E ) = 1.12292 . . ., and the Bessel function of the first kind J 0 (x) with J 0 (0) = 1.
Here and in what follows, the superscript F is attached to process specific quantities in order to distinguish them from universal ones. The symbol ⊗ indicates the convolution in the following sense: 6 Throughout this paper, parameters that are not crucial for the discussion will be suppressed in function arguments. Note that we refrain from including the spin correlation functions G introduced in Ref. [71] here and in what follows, since they are not required at the order considered in this paper.
where f j (x, q) denotes the parton i with momentum fraction x of the proton, and evaluated at momentum transfer q.
The central element of the resummation formula is the so-called Sudakov form factor which resums logarithms of the form L = ln(b 2 M 2 /b 2 0 ), while α s L is treated as being of order unity. The order of the expansion in the exponent then defines the logarithmic accuracy. At leading logarithmic level, only g (1) c has to be taken into account, at NLL also g (2) c and so forth. We give their functional expressions up to the required order in this paper (i.e. g (1) c and g (2) c ) in Appendix A. Clearly, there is a certain amount of freedom in the separation between the "hard" and the "soft" region which can be parametrized by the so-called resummation scale Q. Unless indicated otherwise, we have set Q ≡ M throughout this section; the generalization of the formulas to Q = M and consequently L = ln(Q 2 b 2 /b 2 0 ) can be found in Ref. [43]. In fact, the choice of the resummation scale for the gluon fusion process will be one of the central issues of this paper and will be discussed in more detail in Section 3.
The Born factorσ F,(0) cc in Eq. (5) is given by the parton level cross section at LO. In general, the sum over c accounts for all LO subprocesses that can produce the considered colorless particle. In the gluon fusion process though, only c = g is relevant. An explicit analytical expression forσ F,(0) gg for this process can be found in Eq. (21) of Ref. [68], for example. 8 The resummation coefficients in Eq. (5) can be expanded perturbatively: NLL accuracy requires the knowledge of these coefficients to first order in α s . Evidently, δ(1−z) terms in C cj (equivalently, a particular choice of H F c for one process) defines what is called a resummation scheme [43]. Within a particular resummation scheme, the coefficients C cj can be considered universal, while H F c is process dependent. The entire process dependence in Eq. (5) within a given resummation scheme is thus contained in H F c and σ F,(0) cc . We give the resummation coefficients in the gg → φ scheme (φ ∈ {h, H, A}), which is defined by setting for this process. 8 In the notation of Ref. [68], it is σ The C-coefficients in the gg → φ scheme are [74] C (1) where again q ∈ {u, d, s, c, b} and ζ 2 ≡ ζ(2) = π 2 /6, with Riemann's ζ function. A φ,virt g denotes the finite part of the virtual corrections as defined in Eq. (38) of Ref. [74], i.e.
It can be shown that the unitarity constraint of Eq. (4) can be imposed by replacing in Eq. (7). In addition, this replacement reduces unjustified resummation effects at high transverse momenta, since L vanishes in the limit b → 0 (i.e. p T → ∞), while the large-b limit (small p T ) is preserved.
With the replacement in Eq. (12) the resummed cross section dσ (res) /dp 2 T becomes explicitely Q dependent; however, this dependence formally cancels between [dσ (res) /dp 2 T ] l.a. and [dσ (res) /dp 2 T ] f.o. in Eq. (2). Any residual dependence of the final result is beyond the specific logarithmic order under consideration. The variation of the cross section with Q will be used to estimate the uncertainty due to missing terms of higher logarithmic accuracy.

Components to the matched-resummed cross section
The goal of this paper is to determine p T spectra of neutral MSSM Higgs bosons produced via gluon fusion by matching the NLO result to the resummed NLL approximation.
The relevant NLO matrix elements are taken from Ref. [68], which include the SM-like contributions as well as sbottom, stop and gluino effects (see Fig. 1 for some sample Feynman diagrams). The LO diagrams, e.g. Fig. 1 (2) at p T > 0 9 is governed by the real emission diagrams like the ones shown in Fig. 1 (h) and (i) (and similar ones with quark loops replaced by squark loops). Finally, the virtual diagrams, e.g. Fig. 1 gg , as can be seen from Eq. (10). These contributions allow to calculate [dσ (res) /dp 2 T ] l.a.=NLL .
The expansion of dσ (res) /dp 2 T with respect to α s determines the logarithmic terms at NLO in Eq. (2). The explicit expression can be found in Eq. (72) of Ref. [43], with the corresponding coefficients in Eq. (63) and (64) of that paper.
The resummed expression [dσ (res) /dp 2 T ] NLL has been calculated with a modified version of the program HqT [43][44][45], which determines the NNLO+NNLL p T distribution for the gluon fusion process using the approximation of an infinitely heavy top quark. We modified it for our purposes and implemented the resummation coefficients of Eqs. (10) to include the MSSM effects. Figure 1: A sample of Feynman diagrams for gg → φ contributing to the NLO cross section; (a-c) LO, (d-g) virtual and (h-i) real corrections. The graphical notation for the lines is: solid straight = quark; spiraled = gluon; dashed = scalar (squark or Higgs); spiraled with line = gluino.

Choosing the resummation scale
While the matched cross section is formally independent of the resummation scale Q, the actual numerical result can be quite sensitive to its particular choice due to the truncation at finite logarithmic order. It is therefore vital to determine an "optimal" choice for this scale. The resummation scale Q can be viewed as a scale up to which resummation is extended. The soft and collinear approximation can be trusted only up to a finite value of p T which is determined by a characteristic external scale of the problem. Consequently, there is a maximum value of p T above which resummation is not valid and therefore, Q should not be chosen beyond this value. Due to the constraint of Eq. (4), a too large value of Q not only spoils the prediction for p T < Q, but also affects the large-p T region 10 . The distribution can thus deviate significantly from the fixed-order prediction even in regions where the latter should provide a good approximation.
For the top-quark induced gluon-Higgs coupling, the characteristic scale is of the order of the Higgs mass (m φ ). Consequently, a reasonable range for phenomenological studies is Q ∈ [m φ /4, m φ ], for example. The need for precise predictions requires one to take into account also bottom-quark induced effects to the gluon-Higgs coupling though [49,50]. It has been shown [47] that terms ∼ ln(m b /p T ) appear in the amplitude, which are potentially large and thus could spoil the soft and collinear approximation for p T > m b . This has been used as an argument to choose Q = m b [47], i.e. to effectively turn off Sudakov resummation at p T m b , even though their actual impact has not been studied quantitatively in this context. However, these logarithms are not of Sudakov type as they vanish when p T → 0. In fact, they are closely related to logarithms ln(m b /m h ) which induce an uncertainty already at the level of the total cross section [27]. For a related quantity, namely the cross section with a veto on jets with p jet T > p jet T,veto , it was argued [51] that the impact of the analogous terms ∼ ln(m b /p jet T,veto ) remains moderate, and that one can treat these logarithms as a "finite remainder" together with all other finite terms (power corrections in p jet T,veto ). A similar argument should apply also to the Higgs' p T spectrum.
In this section, we will formulate a pragmatic though quite general way to set the resummation scale. Rather than providing an ex-ante value of Q for a particular process, our method is trial-and-error based and relies on simple expectations on the properties of the matched p T distribution. Roughly speaking, we determine the value of Q as large as possible while requiring that the large-p T behavior of the matched distribution stays reasonably close to the fixed-order prediction.
More precisely, for Higgs masses up to m φ = 300 GeV, we determine Q max as the maximum value of Q for which the resummed p T -distribution stays within the interval [0,2]·[dσ/dp The restriction to the latter interval is needed because, on the one hand, resummation effects are expected to be large for smaller values of p T ; on the other hand, the numerical accuracy of our implementation of the resummation formula becomes unreliable above certain values of p T . The specific value of p max T needs to be chosen case by case. For m φ = 125.6 GeV, it is p max T ≈ 400 GeV, for m φ = 300 GeV, we use p max T ≈ 650 GeV. 10 By "large-pT ", we mean transverse momenta "of the order of the characteristic scale" and beyond.
Neglecting squark effects for the moment, we apply this approach independently to the purely top and bottom induced contributions to the cross section, as well as to the topbottom interference term. Fig. 2 shows these three contributions to the resummed p Tdistribution for different resummation scales in the case of light Higgs production (m h = 125.6 GeV), normalized to the respective fixed-order distribution. The curves for a heavy Higgs of m H = 300 GeV are shown in Fig. 3; those for a pseudo-scalar Higgs of the same mass are very similar to the latter, so we refrain from showing them here.
Larger Higgs masses correspond to a harder p T spectrum since the larger scale of the process leads to less soft gluon radiation. Nevertheless, for m φ = 800 GeV, the numerical accuracy of dσ (res) /dp 2 T becomes unreliable already at p T 700 GeV, so the above procedure for choosing Q cannot be applied. We are therefore forced to modify our criterion for m φ = 800 GeV; our choice is to require |[dσ (res) /dp 2 T ]/[dσ/dp 2 T ] f.o. − 1| = 1/2 at p T = 700 GeV. The corresponding curves for a heavy Higgs of m H = 800 GeV are shown in Fig. 4; again, those for a pseudo-scalar Higgs of the same mass are very similar, so we refrain from showing them here.
Our central scale choice is then defined as Q 0 = Q max /2, while the associated uncertainty is determined by varying Q within the interval [Q 0 /2, 2Q 0 ] (with an additional damping factor for large p T , see Section 4.2). The results of this procedure for a hadronic center of mass energy of √ s = 13 TeV are listed in Tab. 1. Somewhat reassuringly, for m h = 125.6 GeV our value for Q 0,t agrees rather well with the default choice Q = m h /2 of Ref. [45]. On the other hand, our interval for Q 0,b extends to significantly larger values as the one argued for in Ref. [47]. This is even more so for the interference term for which, in our case, the central resummation scale is almost the exact average of the Q 0,t and Q 0,b , while Ref. [47] fully attributed this term to the bottom contribution.
Our result for Q 0,int agrees very well with what was found for the case of jet-veto in Ref. [51]. By analyzing the finite remainder of the bottom contribution, which includes the top-bottom interference in their case, they find Q ≈ 35 GeV to be an appropriate scale  Even though our approach to determine Q seems very pragmatic, the underlying idea is physical, of course. A too large resummation scale Q would overemphasize the Sudakov contribution, typically overshooting the cross section for p T Q. Due to the constraint of Eq. (4), the only way to compensate for this effect is to reduce the cross section at larger transverse momenta, such that it may even become negative. Therefore, by demanding resummation scales that lead to satisfactory matching at high transverse momenta, one indirectly restricts resummation to regions where the soft and collinear factorization is valid. More precisely, we can expect Q to be close to the upper boundary of the range allowed by factorization, certainly not far above that. Note also that there is a certain amount of freedom how this range is defined.
For our later discussion, it will be useful to study the impact of the top-, bottom-and their interference contribution on the shape of the p T distribution, i.e.
1 σ dσ dp T ≡ dσ dp T , with dp T dσ dp T = 1 . In all cases, the bottom-quark distribution is significantly softer than the top contribution, but the difference between the two decreases for larger Higgs masses. This behavior is expected since soft radiation off the quark loop becomes larger for a smaller quark mass, or, equivalently, larger Higgs mass [75].
The shape of top-bottom interference term experiences a number of qualitative and quantitative changes as the Higgs mass increases. For m φ = 125.6 GeV, it can lead to quite some deviations from a pure top-or bottom-dominated shape, see Fig. 5 (a). Whether the spectrum becomes harder or softer depends on the sign of the interference (and thus also on the sign of the Yukawa couplings). 12 Note also that there is a sign change at about p T = 80 GeV. At m φ = 300 GeV, Fig. 5 (b), the qualitative behavior remains roughly the same, but appears to be less distinct. At m φ = 800 GeV, Fig. 5 (c), on the other hand, the shape of the interference term is almost indistinguishable from the top contribution, except for very small p T .
We observe a nice convergence of the resummed and the fixed-order distributions towards large p T , as required by our determination of the matching scale. The curves for a pseudoscalar Higgs are very similar to the scalar case which is why we refrain from showing them here.
As squark effects are typically small due to the fact that squark masses are expected to be of the order of a few hundred GeV, we do not determine separate resummation scales for them. We therefore split the cross section into three terms: • The pure-b contribution is proportional to the square of the bottom-Higgs coupling: σ pure-b ∼ y 2 b . Note that σ pure-b does not include sbottom effects. • The int-b contribution is linearly proportional to the bottom Yukawa coupling: σ int-b ∼ y b . It therefore contains interference terms of the bottom-with the topand squark-loop induced amplitude.
• The no-b contribution is defined as the cross section for y b = 0 and contains topand squark-loop induced terms.
For the pure-b, the int-b, and the no-b contribution, we use the resummation scales Q 0,b , Q 0,int , and Q 0,t of Tab. 1, respectively, despite the fact that these scales were determined by disregarding squark effects. We have checked though that the numerical value of Q 0,int and Q 0,t are hardly affected when squark effects are taken into account.

Input parameters
We present results for the resummed transverse momentum distribution of neutral Higgs bosons produced at the LHC via gluon fusion in various scenarios of the MSSM. Our goal here is not a detailed and comprehensive study of the p T -spectrum in each of these scenarios though. Rather, we will make use of specific scenarios to highlight various features and dependences of the p T -spectrum. Our default choice for the center-of-mass energy is 13 TeV. The central factorization and renormalization scale is set to µ F = µ R = m φ /2. The choice for the central resummation scale is more subtle and is given in Tab. 1. All numbers are obtained with the NLO PDF set of MSTW2008 [76], which implies that the input value for the strong coupling constant is taken as α s (m Z ) = 0.12018. We use the on-shell top and bottom mass with numerical values m t = 173.2 GeV and m b = 4.92 GeV both for the internal propagators and the Yukawa couplings. Terms enhanced by tan β are implicitely resummed [77][78][79][80][81] by reweighing the bottom Yukawa coupling as described in Ref. [82]. Similarly, the stop and the sbottom masses and mixing angles are renormalized as in Ref. [82], in accordance with the definition of the benchmark scenarios to be described in the next section.

MSSM parameter points
We compare results for various MSSM benchmark scenarios, as defined in Ref. [83] 13 and refer to that paper for further details. These benchmark scenarios require the choice of m A scenario m A /GeV tan β m h /GeV m H /GeV  Table 2: Parameter points considered in this paper. The full definition of the scenarios is given in Ref. [83]; for the light-stop scenario, however, we modify the soft SUSY breaking wino and bino mass terms as well as the µ-parameter as suggested in Ref. [82], where m 2 = µ = 400 GeV and m 1 = 340 GeV, in order to evade constraints on the stop and sbottom masses presented by ATLAS and CMS. The particular parameter points defined here will be refered to in the text as "scenario(m A /GeV,tan β)"; for example, the first parameter point in the table is τ -phobic(800,16) in this notation. The Higgs masses are evaluated with FeynHiggs [87-95] which we also apply to determine the the corresponding Higgs couplings in the various scenarios. and tan β. Using the exclusion plots of Ref. [83] and HiggsBounds [84][85][86], we identified proper (i.e. not yet excluded) parameter choices within the m A -tan β plane, while requiring that m h = 125.6 ± 0.7 GeV (except for the light-stop scenario). The scenarios used for our analysis in Section 5 are defined in Tab. 2.

Theoretical uncertainties
The main sources of theoretical uncertainty on our result for the p T distribution are due to missing higher order effects, as well as the uncertainty from the PDFs and α s (m Z ). The latter two are usually estimated by following the so-called PDF4LHC recipe [96]. They will, however, not be part of our analysis within this manuscript.
The former are typically estimated by a variation of unphysical scales that emerge at  finite perturbative or logarithmic order. In our case, these are the renormalization, the factorization, and the resummation scale. While for µ F and µ R , we follow the standard procedure of considering the maximum variation of the cross section when 2µ R /m φ and 2µ F /m φ are taken from the set {1/2, 1, 2}, while excluding the values for which µ R /µ F ∈ {1/4, 4}. The impact of the choice of the resummation scale, on the other hand, we estimate by varying Q/Q 0 within the interval [1/2, 2]. However, a variation within this region at large p T would grossly overestimate the uncertainty at p T m φ , where the prediction should be well described by the fixed-order distribution. We therefore modulate the error band resulting from Q-variation by a damping factor which effectively switches off the Q-uncertainty for p T m φ . Finally, we add the uncertainty estimated from µ F -and µ R -variation and the one induced by Q in quadrature.

Results
We are now ready to present our results for the transverse momentum spectrum of MSSM Higgs bosons produced in gluon fusion through NLO+NLL. b For reference, Fig. 6 shows the p T -spectrum of a SM Higgs boson of mass m h = 125.6 GeV. (dσ x /dp T ) / (dσ htl /dp T ) p T [GeV] pp @ 13 TeV b+t t Figure 7: Exact top and bottom mass dependence of the transverse momentum distribution at NLO+NLL. In the dashed curve, bottom quark effects are set to zero. The normalization of these curves is to the respective LO total cross section times the NLO+NLL result in the heavy top limit.
The resummed cross section is finite in the limit p T → 0 and smoothly matches the fixed order curve at large transverse momenta (p T m h ). The uncertainty band is obtained through scale variation of µ F , µ R and Q, following the procedure described in Section 4.2. In order to compare to other calculations, it may be useful to consider the ratio of the p T distribution which includes the full quark mass dependence to the result in the heavy-top limit (reweighted by the full LO inclusive cross section for gg → H). The corresponding curve for the SM is shown in Fig. 7 and can be compared to analoguous plots of Refs. [47,49,50,97]. Disregarding the specific normalization in these papers, the behaviour of the curve which includes both top-and bottom-quark effects is quite different in the various approaches, in particular towards small values of p T . For example, in Ref. [50], where a common resummation scale for the top-and the bottom-quark effects of Q 0,t = Q 0,b = Q 0,int = m φ /2 was chosen, the curve drops only by about 6% between p T = 100 GeV and p T = 0. With a separate resummation scale for the bottom-effects of Q 0,b = Q 0,int ∈ [m b , 4m b ] as suggested in Ref. [47], this effect becomes much more pronounced and amounts to about (27 ± 9)%. In the POWHEG approach of Ref. [49], on the other hand, the drop in the curve is roughly 20%, while the MC@NLO [98] result of Ref. [97] with a drop of 5% is quite similar to the analytic resummation.
Separate resummation scales for the top, the bottom, and the interference term as given in Tab. 1, on the other hand, lead to a drop of 11%, which is of a size as expected considering the magnitude of the new scale choices compared to the ones of the previous studies in Refs. [47,50].
It has been shown that differences in the various approaches (analytic resummation, POWHEG, MC@NLO) become much smaller by a simultaneous adjustment of the corresponding intrinsic scales (Q, h fact , shower scale). Due to the similarities in their NLO matching, this leads to an excellent agreement for various scale choices [99] for the analytic resummation and MC@NLO. Also, the initially observed large differences to the POWHEG approach are alleviated [100] at least when in all approaches the scale for the bottom contribution is choosen of the order of the bottom mass.
The p T distributions of the light Higgs boson in the various scenarios of Tab. 2 are virtually indistinguishable from the SM distribution shown in Fig. 6. This is because the observation of a Higgs particle at about m h = 125 GeV typically constrains the parameter space of the MSSM in such a way that the light Higgs is SM-like.
In order to quantify the deviations between the MSSM and the SM prediction, Fig. 8 (a) shows the ratios of the resummed p T distributions at NLO+NLL for the various MSSM scenarios S with respect to the SM one. The difference to the SM is typically at the 1-3% level; only scenario m max h (300,6.5) deviates by up to 9%. All curves are below one, because their respective total cross sections are smaller than the SM one. Considering the ratio for the shapes, i.e. the normalized distributions, Eq. (13), we find variations at the 2%-level, see Fig. 8 (b), with the only exception again m max h (300,6.5) which, however, still stays within 4% of the SM prediction. Apparently, the harder spectrum for the MSSM scenarios compared to the SM is due to a slightly larger negative int-b term, recall Fig. 5. The numerical effects observed here are roughly of the same size as those observed in Ref. [49] (see the right plot of Fig. 8 in that paper). The fact that for all scenarios S, the ratio N S = 1 occurs at roughly the same value of p T ≈ 30 GeV is a consequence of the similar "barycenter" p T = σ −1 tot dp T p T (dσ/dp T ) of the distributions. While the predictions for the light Higgs are very SM-like, this is not the case for the heavy and pseudo-scalar Higgs, see Fig. 9. 14 We show curves for various scenarios with m A = 800 GeV, where m H ≈ m A . Clearly, the absolute size of the cross section for both H and A depends strongly on the respective scenario and the value of tan β (see also Ref. [82]). Indeed for the heavy Higgs, in all scenarios the cross section increases with the value of tan β, which is caused by the fact that the bottom contribution strongly increases and eventually becomes the dominant contribution to the cross section. One remarkable observation is that for the pseudo-scalar Higgs the curves in the τ -phobic scenarios for the two different values of tan β are quite close, in contrast to all other scenarios. In fact, at large values of p T , the cross section for tan β = 16 is even bigger than the one for tan β = 29.5. The reason for this behavior will be discussed further below.
In Section 3, we split the cross section into the three contributions pure-b, no-b, and intb for which separate resummation scales were determined. The relative contribution of these three terms to the heavy Higgs p T distribution for (a) the τ -phobic(800,16) and (b) the τ -phobic(800,29.5) scenario is shown in Fig. 10. Note that by definition the pure-b (red, solid), no-b (blue, dashed) and int-b curve (brown, dash-dotted) add up to one (black, dash-double dotted). For comparison, we also include a curve for the "pure-t contribution" (green, dotted) which is defined to be proportional to the square of the top-Higgs coupling y t .
For τ -phobic(800,16), we find a rather large cancellation between the positive no-b-and pure-b-, and the negative int-b term, see Fig. 10 (a). It shows the importance of the proper treatment of the int-b term in the resummation procedure and justifies a separate resummation scale as introduced in Section 3. By comparing to the pure-t contribution, we also observe that the squark effects are of the order of the overall contribution and therefore very relevant.
Going to tan β = 29.5, the cancellations among the individual contributions are less severe,  (dσ x /dp T ) / (dσ sum /dp T ) p T [GeV] tauphobic, m A = 800 GeV, pp @ 13 TeV . The pure bottom term (red, solid), the no-b term (blue, dashed) and their interference (brown, dash-dotted) add up to one, which is marked for reference (black, dash-double dotted). For comparison, the pure top contribution is shown as well (green, dotted). see Fig. 9 (b). As expected, the cross section is largely dominated by bottom-quark effects, i.e., the pure-b and the int-b term. These results substantiate that the contribution of the bottom loop causes the increase of the cross section that we observed in Fig. 9 (a) at high tan β.
Let us now compare these observations in the τ -phobic to the m mod+ h scenario shown in Fig. 10   ). The structure of the various contributions to the cross section is quite different from the one for the heavy Higgs. Since the int-b term is positive here, all curves remain between 0 and 1. In both scenarios, for the smaller value of tan β, each term contributes at least about 20% to the cross section, and none of them exceeds 45%. The largest contribution is due to the int-b term for most transverse momenta. At high tan β, the pure-b contribution becomes clearly dominant again. While the int-b term remains sizable, both pure-t and no-b terms are negligible, especially in the m mod+ h scenario.
With the results of Fig. 11 (a) and (b), it is interesting to take another look at the behavior of the distribution for the pseudo-scalar Higgs shown in Fig. 9 (b) in the τ -phobic scenarios. The splitting into the individual contributions suggests that there is no deep reason for the similarity of the curves for tan β = 16 and tan β = 29.5. The hierarchy of the various contributions to the cross section in the τ -phobic scenarios is not very different from the m mod+ h scenarios of Fig. 11. It rather seems to be an accidental interplay of the top-and bottom-quark effects so that the absolute size of the increase of the pure-b contribution from tan β = 16 to tan β = 29.5 is compensated by a decrease of similar size of the no-b and int-b term.
To finalize the analysis of the p T distributions shown in Fig. 9, we study their shapes in the various scenarios for the heavy Higgs in Fig. 12 (a) and for the pseudo-scalar Higgs in Fig. 12 (b) by considering again the ratio of shapes defined in Eq. (16). Note that the normalization is for a "SM Higgs" of mass 800 GeV. For the heavy Higgs in Fig. 12 (a), we observe generally small deviations between the curves. While the biggest difference occurs at large transverse momenta, the similarity in shape at small p T is remarkable. Their deviation for the SM curve is quite large though, reaching up to 60%, and clearly showing the dominance of the pure-b term by the significantly softer spectrum, see Fig. 5. For most scenarios, however, the softness decreases with increasing tan β, with the exception of the τ -phobic scenarios. This is again the impact of the negative int-b term. The softening of the spectrum due to an enhanced b-contribution is in agreement with the observations of (dσ x /dp T ) / (dσ sum /dp T ) p T [GeV] tauphobic, m A = 800 GeV, pp @ 13 TeV (dσ x /dp T ) / (dσ sum /dp T ) p T [GeV] mhmodp, m A = 800 GeV, pp @ 13 TeV Figure 11: Same as Fig. 10, but for the pseudo-scalar Higgs.  Refs. [49,75].
Considering the pseudo-scalar Higgs in Fig. 12 (b), the spread of the curves is significantly larger both at small and high transverse momenta, leading to a more enhanced difference in shape of the resummed p T distributions in the various scenarios. Since the interference contributions are strictly positive in this case, the corrections for all scenarios, including τ -phobic, decrease with increasing tan β. We also note that the m mod+ h and m mod− h curves are practically indistinguishable in this case.
As a final study, we compare the transverse momentum distributions at different Higgs masses. Fig. 13 shows the ratio of the shapes for the heavy scalar in the m max h (300,6.5) and in the m max h (800,6.5) scenario. For comparison, the same ratio is shown at fixed order (blue, dashed). In addition, the ratio of the p T -shape of a SM Higgs at 125.6 GeV and the heavy Higgs in the m max h (800,6.5) scenario (green, dotted) is given. The spectra at low Higgs masses are significantly softer due to increased soft radiation. This observation is consistent with the behavior of Figs. 6 and 9 (a), which differ by an order of magnitude in the difference between the minimum at p T = 300 GeV and the maximum of the curve. Furthermore, Fig. 13 confirms that the shape at high transverse momenta is driven by the fixed order cross section, as expected. The harder spectrum at high Higgs masses is caused by the fact that the colliding gluons carry more energy in a production of a heavy particle, which makes it more likely to emit harder gluons and therefore, to produce a harder Higgs boson.

Conclusions
We have presented typical MSSM effects on the p T spectrum of a neutral Higgs boson produced at the LHC. Special emphasis has been put on the impact of bottom quarks. While the current experimental data imply a small, SM-like bottom-Yukawa coupling for the light Higgs boson, the production mechanism for the heavy and the pseudo-scalar Higgs can be dominated by bottom-quark loops.
Through a pragmatic argumentation based on simple theoretical and physical expectations, we derived separate resummation scales for the pure-b and int-b contributions which turn out to be significantly larger than the bottom-quark mass, and smaller than the Higgs boson mass.
We find the well-known behavior that the bottom loop typically softens the p T -spectrum. For the pseudo-scalar Higgs boson, the int-b term is typically positive and leads to a further softening of the spectrum. For the light and heavy CP-even Higgs bosons, on the other hand, the int-b term is typically negative, and makes the spectrum a bit harder. The latter effect is specifically visible for the light Higgs where the spectrum in all MSSM benchmarks is harder than in the SM, because the relative contribution of the negative int-b term is larger. In all scenarios, we observe an enhanced importance of the bottomquark contributions for the heavy and pseudo-scalar Higgs, which become by far dominant at large values of tan β. Indeed, the corresponding spectra are clearly softer than in the SM. Finally, we confirmed that larger Higgs masses lead to reduced soft gluon emission and therefore a harder spectrum.
The resummed p T -distributions through NLO+NLL have been implemented in the program MoRe-SusHi, which advances the program SusHi to small-p T distributions. The code is publically available and can be found on the SusHi homepage. 15 As an outlook, one may consider combining the consistent NLO+NLL results presented in this paper for the p T -distribution of MSSM Higgs bosons with the NNLO+NNLL distribution in the infinite-top mass approximation, and possibly even with the p T spectrum of the Higgs boson produced through bottom-quark annihilation. This is certainly beyond the scope of this paper and is left for a future publication. account. Their functional expressions for Q = µ R = M are the following: where and β 0 = (11 C A − 2 N f )/12 and β 1 = (17 C 2 A − 5 C A N f − 3 C F N f )/24 are the first two coefficients of the β function.
The process and resummation scheme independent coefficients needed at NLL have been known for some time for gluon-induced processes [40]; they read where C F = 4/3, C A = 3, and N f = 5 is the number of active quark flavors.