Higgs production via gluon fusion in the POWHEG approach in the SM and in the MSSM

We consider the gluon fusion production cross section of a scalar Higgs boson at NLO QCD in the SM and in the MSSM. We implement the calculation in the POWHEG approach, and match the NLO-QCD results with the PYTHIA and HERWIG QCD parton showers. We discuss a few representative scenarios in the SM and MSSM parameter spaces, with emphasis on the fermion and squark mass effects on the Higgs boson distributions.


Introduction
Understanding the mechanism that leads to the breaking of the electroweak symmetry and that is responsible for the generation of the mass of the elementary particles is one of the major challenges of high energy physics. The search for the Higgs boson(s) is currently under way at the Tevatron and at the LHC, and limits on the Higgs mass spectrum have already been set [1,2]. This search requires an accurate control of all the Higgs production and decay mechanisms, including the effects due to radiative corrections [3].
In the Standard Model (SM) the gluon fusion process [4] is the dominant Higgs production mechanism both at the Tevatron and at the LHC. The total cross section receives very large next-to-leading order (NLO) QCD corrections, which were first computed in ref. [5] in the so-called heavy-quark effective theory (HQET), i.e. including only the top-quark contributions in the limit m t → ∞. Later calculations [6,7,8,9,10] retained the exact dependence on the masses of the top and bottom quarks running in the loops. The next-to-next-to-leading order (NNLO) QCD corrections are also large, and have been computed in the HQET in ref. [11]. The finite-top-mass effects at NNLO QCD have been studied in ref. [12] and found to be small. The resummation to all orders of soft gluon radiation has been studied in refs. [13,14]. Leading third-order (NNNLO) QCD terms have been discussed in ref. [15]. The role of electroweak (EW) corrections has been discussed in refs. [16,17,18,19]. The impact of mixed QCD-EW corrections has been discussed in ref. [20]. The residual uncertainty on the total cross section depends mainly on the uncomputed higher-order QCD effects and on the uncertainties that affect the parton distribution functions (PDF) of the proton [21,22,3].
The Higgs sector of the Minimal Supersymmetric Standard Model (MSSM) consists of two SU (2) doublets, H 1 and H 2 , whose relative contribution to electroweak symmetry breaking is determined by the ratio of vacuum expectation values of their neutral components, tan β ≡ v 2 /v 1 . The spectrum of physical Higgs bosons is richer than in the SM, consisting of two neutral CP-even bosons, h and H, one neutral CP-odd boson, A, and two charged bosons, H ± . The couplings of the MSSM Higgs bosons to matter fermions differ from those of the SM Higgs, and they can be considerably enhanced (or suppressed) depending on tan β. As in the SM, gluon fusion is one of the most important production mechanisms for the neutral Higgs bosons, whose couplings to the gluons are mediated by top and bottom quarks and their supersymmetric partners, the stop and sbottom squarks.
In the MSSM, the cross section for Higgs boson production in gluon fusion is currently known at the NLO. The contributions arising from diagrams with quarks and gluons can be obtained from the corresponding SM results [6,7,8,9,10] with an appropriate rescaling of the Higgs-quark couplings. The contributions arising from diagrams with squarks and gluons were first computed under the approximation of vanishing Higgs mass in ref. [23], and the full Higgs-mass dependence was included in later calculations [8,9,10,24]. The contributions of two-loop diagrams involving top, stop and gluino to both scalar and pseudoscalar Higgs production were computed under the approximation of vanishing Higgs mass in refs. [25,26], whose results were later confirmed and cast in a compact analytic form in refs. [27,28].
The approximation of vanishing Higgs mass can provide reasonably accurate results as long as the Higgs mass is well below the threshold for creation of the massive particles running in the loops. For the production of the lightest scalar Higgs, this condition does apply to the two-loop diagrams involving top, stop and gluino, but it obviously does not apply to the corresponding diagrams involving the bottom quark, whose contribution can be relevant for large values of tan β. In turn, the masses of the heaviest scalar and of the pseudoscalar might very well approach (or exceed) the threshold for creation of top quarks or even of squarks. Unfortunately, retaining the full dependence on the Higgs mass in the quark-squark-gluino contributions has proved a rather daunting task. A calculation based on a combination of analytic and numerical methods was presented in ref. [29] (see also ref. [30]), but neither explicit analytic results nor a public computer code have been made available so far. However, ref. [31] presented an approximate evaluation of the bottom-sbottom-gluino contributions to scalar production, based on an asymptotic expansion in the large supersymmetric masses that is valid up to and including terms of where m φ denotes a Higgs boson mass and M SU SY denotes a generic superparticle mass. An independent calculation of the bottom-sbottom-gluino contributions, restricted to the limit of a degenerate superparticle mass spectrum, was also presented in ref. [32], confirming the results of ref. [31]. More recently, ref. [28] presented an evaluation of the quark-squark-gluino contributions to pseudoscalar production that is also based on an asymptotic expansion in the large supersymmetric masses, but does not assume any hierarchy between the pseudoscalar mass and the quark mass, thus covering both the top-stop-gluino and bottom-sbottom gluino cases.
The total cross section, without acceptance cuts, provides important information about the Higgs boson production rate. On the other hand, especially at the LHC, it is very likely that a Higgs boson is produced in association with a jet, generating a transverse momentum p H T of the Higgs boson. The first studies on Higgs+jet final states in the SM were performed in ref. [33], using the real-parton emission amplitudes that enter the calculation of the NLO-QCD corrections to the inclusive Higgs production cross section. The NLO corrections to the Higgs+jet final state at large transverse momentum of the Higgs boson were subsequently studied in ref. [34], and the resummation of all the logarithmically enhanced terms, matched with the NLO calculation of the Higgs+jet final state, was discussed in ref. [35]. The impact on the Higgs+jet final state of the NLO-EW corrections and of the finite masses of the particles running in the loops was discussed in ref. [36]. In the MSSM, the production of a neutral Higgs boson in association with one jet was discussed in refs. [37,38].
A different set of observables are the differential distributions of the Higgs boson, inclusive over QCD radiation. For the SM case, results were presented in ref. [39] at NLO-QCD and in ref. [40] at NNLO-QCD. The Higgs boson transverse momentum spectrum, including the NLO-QCD corrections matched with the resummation of next-to-next-to-leading logarithmic (NNLL) enhanced terms, has been studied in refs. [41,42,43]. In the MSSM a study of the Higgs distributions, at NLO-QCD, was discussed in ref. [32]. If a new scalar particle is discovered at the Tevatron or at the LHC, a major question will be to determine whether it is a Higgs boson and, in that case, whether it belongs to the particle spectrum of the SM, of the MSSM or of any other model. An example could be represented by a MSSM Higgs boson whose production cross section is close to the production cross section for a SM Higgs boson of equal mass. In this case, an accurate study of the differential distributions involving the Higgs boson might shed some light on the underlying model. A precise analysis of the experimental data requires the use of NLO-QCD results merged with the description of initial-state multiple gluon emission via a QCD parton shower (PS) Monte Carlo (MC) such as HERWIG [44] or PYTHIA [45]. However, the merging of NLO-QCD matrix elements with PS faces the problem of avoiding double counting, as addressed in refs. [46,47]. The POWHEG method [48] allows to systematically merge NLO calculations with vetoed PS, avoiding double counting and preserving the NLO accuracy of the calculation. The procedure can be implemented using a set of tools and results available in the so-called POWHEG BOX [49]. The latter provides a general framework that exploits the universal nature of initial-state collinear divergences and the factorization property of soft radiation to automatize the subtraction of all the soft and/or collinear divergent terms from the NLO matrix elements of an arbitrary process. The POWHEG method does not rely on the details of the shower MC and, by construction, guarantees an accuracy at NLO + leading logarithmic (LL) QCD. In ref. [51] it has been shown that, with an appropriate choice of scale for the strong coupling constant, the merging procedure can also reproduce the next-to-leading logarithmic (NLL) terms.
At present, no code exists that merges the NLO-QCD results for the gluon fusion process with a QCD PS, retaining the exact dependence on the Higgs mass and on the masses of the particles running in the loops. Two implementations of the NLO-QCD results merged with a PS are available 1 : the one in MC@NLO [50] and the one in the POWHEG BOX framework [51]. However, both implementations are limited to the SM case, and beyond LO they only include results computed in the HQET. Conversely, the codes HIGLU [53] and iHixs [54] contain the full dependence on the Higgs and quark masses up to NLO, and HIGLU also allows to include the full squark-gluon contributions from ref. [24], but neither code is matched to a shower MC. Recently, a step toward the inclusion of the finite-quark-mass effects in PS was taken in ref. [55], where parton-level events for Higgs production accompanied by zero, one or two partons are generated with matrix elements computed in the HQET, and then, before being passed to the PS, they are re-weighted by the ratio of the exact one-loop amplitudes over the approximate ones. This procedure is equivalent to generating events directly with the exact one-loop amplitudes, yet it is much faster.
We aim to provide a code that fills the remaining gap, using matrix elements that include the dependence on the masses, both in the SM and in the MSSM, properly matched to an external shower MC. For the SM case, we use NLO matrix elements with full dependence on the Higgs, top and bottom masses. For the MSSM case, we use matrix elements with exact dependence on the quark, squark and Higgs masses in the contributions of real-parton emission. For the two-loop virtual contributions, the approximation of vanishing Higgs mass is employed in the diagrams involving superpartners, while the rest is computed exactly.
The plan of the paper is as follows: in section 2 we describe the basic features of the POWHEG implementation of gg → φ; in section 3 we discuss our SM implementation with exact dependence on the fermion masses, presenting a numerical analysis valid for an on-shell Higgs; section 4 is devoted to analyzing the MSSM case; finally, in section 5 we draw our conclusions.
2 POWHEG implementation of gg → φ In this section we briefly discuss the implementation of the gluon-fusion Higgs production process in the POWHEG BOX framework, following closely ref. [51] (see also ref. [56]). We fix the notation keeping the discussion at a general level, without referring to a specific model. In the next sections the formulae presented below will be specialized to the SM and MSSM cases.
The generation of the hardest emission is done in POWHEG according to the following formula: In the equation above the variablesΦ 1 ≡ (M 2 , Y ) denote the invariant mass squared and the rapidity of the Higgs boson, which describe the kinematics of the Born (i.e., lowest-order) process gg → φ.
The variables Φ rad ≡ (ξ, y, φ) describe the kinematics of the additional final-state parton in the real emission processes. In particular, denoting by k ′ 2 the momentum of the final-state parton in the partonic center-of-mass frame, or we have where s is the partonic center-of-mass energy squared. The factorB(Φ 1 ) in eq. (1) is related to the total cross section computed at NLO in QCD. It contains the value of the differential cross section, for a given configuration of the Born final state variables, integrated over the radiation variables. The integral of this quantity on dΦ 1 without acceptance cuts yields the total cross section. This factor is responsible for the correct NLO-QCD normalization of the result, and is computed in the initialization phase using the real and virtual NLO-QCD corrections.
The terms within curly brackets in eq. (1) describe the real emission spectrum of an extra parton: the first term is the probability of not emitting any parton with transverse momentum larger than a cutoff p min T , while the second term is the probability of not emitting any parton with transverse momentum larger than a given value p T times the probability of emitting a parton with transverse momentum equal to p T . The sum of the two terms fully describes the probability of having either zero or one additional parton in the final state. The probability of non-emission of a parton with transverse momentum k T larger than p T is obtained using the POWHEG Sudakov form factor Finally, the last term in eq. (1) describes the effect of the qq → φg channel, which has been kept apart in the generation of the first hard emission because it does not factorize into the Born cross section times an emission factor.
We now discuss the various terms appearing in eq. (1) in more detail. We have: where with L ab the luminosity for the partons a and b. In eq. (5) " c. r." denotes the collinear remnants multiplied by the relevant parton luminosity. The remnants are the finite leftovers after the subtraction of the initial-state collinear singularities into the parton distribution function is performed, and their explicit expressions are given in eqs. (2.36), (2.37) and (3.7)-(3.10) of ref. [51]. The function B gg (Φ 1 ) in eq. (6) represents the squared matrix element of the Born contribution to the process, averaged over colors and helicities of the incoming gluons, and multiplied by the flux factor 1/(2M 2 ). It is given by where H is the form factor for the coupling of the Higgs boson with two gluons, whose explicit form depends on the particle content of the model considered and will be detailed in the following sections. It is decomposed in one-and two loop parts as The regularized two-loop virtual contributions are contained in In the equation above, µ R and µ F are the renormalization and factorization scale, respectively, C A = N c (N c being the number of colors), and β 0 = (11 C A −2 N f )/6 (N f being the number of active flavors) is the one-loop beta function of the strong coupling. The hatted functionsR ij in eqs. (8)-(10) are the Frixione, Kunst and Signer [57,58] infraredsubtracted counterparts of R iĵ where R ij are the squared amplitudes, averaged over the incoming helicities and colors and multiplied by the flux factor 1/(2s), for the NLO partonic subprocesses (gg → φg, gq → φq, qg → qφ): The complete real matrix elements that enter the POWHEG Sudakov form factor, eq. (4), read where the functions R ab are the non-infrared-subtracted counterparts of eqs. (8)-(10). The probability for the emission of the first and hardest parton is described with the exact matrix element in all the phase space regions. Finally, the contribution of the qq → φg channel is with The functions A gg , A qg in eqs. (15)- (17) and A qq in eq. (21) depend on the particle content of the model considered, and will be defined in the following sections.

SM results
The current public release of POWHEG [49] contains matrix elements evaluated in the HQET. It also gives the user the possibility of rescaling the termB(Φ 1 ) in eq. (1) by a normalization factor defined as the ratio between the exact Born contribution where the full dependence from the top and bottom masses is kept into account and the Born contribution evaluated in the HQET. In the following we describe the modifications we have introduced in the code to include the full fermion-mass dependence at the NLO and the effect of the two-loop EW corrections.

Modifications in POWHEG
The inclusion of the fermion-mass effects is achieved using for the functions H 1ℓ , H 2ℓ , A gg , A qg , A qq the exact results instead of those computed in the HQET. For the Born term we have where T F = 1/2 is the matrix normalization factor of the fundamental representation of SU (N c ), and λ q is a normalization factor for the Higgs-quark coupling. In the SM case λ q = 1 for both the top and the bottom quark.
The form factor H 2ℓ contains the mass-dependent contribution of the two-loop virtual corrections, and can be cast in the following form: where given in terms of harmonic polylogarithms can be found in ref. [9]. It should be noticed that G (2ℓ,C R ) 1/2 depends on the choice of renormalization scheme for the quark mass entering the one-loop part of the form factor. In ref. [9] expressions for G (2ℓ,C R ) 1/2 with on-shell (OS) or MS parameters are presented. In our implementation we allow the choice among the OS, MS or DR renormalization schemes.
Concerning the real emission contributions, we have for the gg → Hg channel where the functions A 2 and A 4 can be cast in the following form: with Explicit expressions for the functions b 1/2 (s q , t q , u q ) and c 1/2 (s q , t q , u q ) are given in ref. [10].
The function A qq (s, t, u) relevant for the qq → Hg channel is given by and d 1/2 (s q , t q , u q ) can be found in ref. [10]. Finally A qg (s, t, u) relevant for the qg → Hg channel can be obtained from The two-loop EW corrections are included as a factor (1 + δ EW ) which multiplies the termB(Φ 1 ) in the first line of eq. (1). This choice follows from the current structure of POWHEG where the qq → Hg channel is kept apart, because it is not proportional to the Born cross section in the collinear limit. In the SM case, the values of the correction δ EW as a function of the Higgs boson mass can be obtained from ref. [18].

SM: numerical results
In this section we present numerical results for the production of an on-shell Higgs boson in the SM. We focus our analysis on the inclusion of the exact quark-mass dependence in the NLO corrections and on the effect of the EW corrections. We also consider the effect of merging POWHEG with a PS. The results have been obtained for the LHC with center-of-mass energy of 7 TeV, using the following numerical values for the physical input parameters: G µ = 1.16637 · 10 −5 GeV −2 , m t = 172.5 GeV and m b = 4.75 GeV [3]. We have used the MSTW2008 [59] NLO set of PDF to describe to partonic content of the proton. In the code the value of α s (m Z ) is set accordingly to the choice made in the PDF set: in our case α s (m Z ) = 0.12018. When discussing the distributions in the Higgs transverse momentum p H T , a cut p H T > 0.8 GeV has been enforced. The renormalization and the factorization scales have been set equal to the Higgs boson mass: In the left panel of figure 1 we plot the total Higgs production cross section, without acceptance cuts, in three different approximations: the dot-dashed line corresponds to the current public POWHEG implementation, in which the NLO-QCD corrections are computed in the HQET and are then rescaled with the exact Born cross section (which includes the full dependence on the top and bottom masses); the dashed line corresponds to the POWHEG implementation presented in this paper, where the complete NLO-QCD calculation is employed, i.e. the top and bottom contributions are treated exactly in the NLO corrections; the solid line also includes the effect of the EW corrections.
In the right panel of figure 1 we plot the ratio between the Higgs production cross section obtained using our version of POWHEG and the one obtained using the current public version. exact treatment of the quark masses results in an increase up to ∼ 6% in the cross section, with a further increase (up to a combined ∼ 10%) when the EW corrections are included. This effect is mainly due to the bottom-quark contribution, which is not negligible when the Higgs boson is light.
For m H 180 GeV the quark-mass effects and the EW corrections have opposite sign, resulting in a ∼ −2% correction in the Higgs mass range up to 2 m t . Above the threshold for real top-quark production, where the approximation m t → ∞ is not valid, the large corrections due to the quarkmass effects in the QCD contribution are partially screened by the EW corrections. In the rest of the section we discuss the kinematic distributions of a SM Higgs boson at NLO QCD. Since the effect of the EW corrections is very close to an overall rescaling of the total cross section, we neglect them in the following and focus on the effect of the QCD corrections.
In figure 2 we show the rapidity distribution for a Higgs boson with mass m H = 120 GeV. In the left panels, the dashed red lines correspond to the distributions obtained using the current implementation of POWHEG, and the solid blue lines correspond to our implementation. In the right panels we plot the ratio between the distributions obtained with the two implementations. We compare two different approximations: in the upper panels we show the pure (i.e., fixed-order) NLO-QCD calculation, while in the lower panels we show the event distributions according to the basic formula of POWHEG, eq. (1), which includes the effects of the Sudakov form factors, merged also with the PYTHIA QCD PS. As appears from the plots, the exact treatment of the quark masses results in a ∼ 5% enhancement, uniformly distributed over the whole rapidity range. In the absence of acceptance cuts, the results obtained by combining POWHEG and PYTHIA do not differ significantly from the pure NLO-QCD results. In figure 3 we show the transverse momentum distribution for a Higgs boson of mass m H = 120 GeV. In the left panels we compare the current POWHEG implementation with ours but, differently from figure 2, we show separately the pure NLO-QCD, POWHEG and POWHEG + PYTHIA calculations. In the right panels we plot the full results (blue, solid lines) and the ones obtained by introducing in POWHEG only the exact top-mass dependence (black, dashed lines), both normalized to the results of the current POWHEG implementation.
The pure NLO-QCD calculation (upper panels) diverges for vanishing Higgs transverse momentum, although in figure 3 the divergence is masked by the lower cut on p H T . The plot on the right shows that for p H T m t it is the inclusion of the bottom-quark contribution in the NLO corrections that gives rise to a positive correction with respect to the result obtained with the current POWHEG version, while for p H T m t there is a substantial modification of the distribution, with a large negative correction driven by the use of the exact top-quark contribution.
The inclusion of the effect of the Sudakov form factor (central panels) modifies the distribution in such a way that the latter vanishes in the limit p H T → 0. As can be read from eq. (1), the probability of emitting the Higgs in association with a parton depends on the product ∆ × R/B, where R is the squared matrix element for real-parton emission, B is the Born amplitude, and ∆ is the Sudakov factor, which in turn is exponentially suppressed by R/B, see eq. (4). In the current POWHEG implementation, the emission probability is computed in terms of the ratio R(t, ∞)/B(t, ∞), where both R and B are evaluated in the HQET. 2 For small p H T , the Sudakov factor with exact top and bottom mass dependence, ∆(t+b, exact), used in our implementation is smaller than the corresponding factor ∆(t, ∞) used in the current POWHEG implementation, because R(t + b, exact)/B(t + b, exact) > R(t, ∞)/B(t, ∞). This inequality holds for two reasons: first, the p H T distribution is proportional to R, and R(t + b, exact) > R(t, ∞) for p H T < 200 GeV [36]; second, the inclusion of the bottom contribution reduces the LO cross section with respect to the result obtained in the HQET [6]. Thus, as shown in the right plot, for small p H T the Sudakov factor suppresses the p H T distribution by almost 10% with respect to the result obtained in the current POWHEG implementation. Since the emission probability is also directly proportional to the ratio R/B, starting from p H T ≃ 30 GeV this factor prevails over the Sudakov factor, and the distribution with exact dependence on the quark masses becomes larger than the one in the current POWHEG implementation by up to ∼ 15%. Finally, for p H T m t the inclusion of the full top-mass dependence leads to a negative correction, similar to the one already observed in the pure NLO-QCD calculation. The inclusion of multiple gluon emission with the PYTHIA QCD-PS (lower panels) does not change dramatically -in the absence of acceptance cuts -the results obtained including only the hardest emission. Figure 4 shows the transverse momentum distribution for a Higgs boson of mass m H = 500 GeV. The meaning of the different curves is the same as in figure 3. In the pure NLO-QCD calculation (upper panels), the effect of using a finite top mass in the NLO corrections is always negative and monotonically decreasing for increasing transverse momentum. Contrary to the light-Higgs case, the use of a finite top mass in the Sudakov form factor (central panels) yields an enhancement of the transverse momentum distribution at small values of p H T . For p H T 80 GeV the mass effects induce instead a dampening of the spectrum. As for the case of m H = 120 GeV, the introduction of the QCD PS (lower panels) does not modify substantially the distributions. As can be seen from the plots on the right, for such a large value of the Higgs mass the bottom contribution is highly suppressed.
In order to illustrate the uncertainty arising from the use of different PS algorithms, in the left panel of figure 5 we compare the transverse momentum distributions obtained with our implementation of POWHEG matched with either PYTHIA (red dashed line) or HERWIG (blue solid line), for a Higgs boson of mass m H = 120 GeV. In the right panel of figure 5 we plot the ratio of the distribution obtained with POWHEG+HERWIG over the one obtained with POWHEG+PYTHIA.     yields a differential cross section larger by up to 10-20% with respect to PYTHIA in the region of small transverse momenta (p H T ≤ 20 GeV), and smaller by 5-10% for large values of p H T . We remark that this behavior is not specific to our implementation of POWHEG. Indeed, the current implementation yields differences of similar size when matched with HERWIG instead of PYTHIA.
In our SM analysis the top and bottom masses are renormalized in the on-shell scheme. Therefore, the dependence of the differential cross section on the renormalization and factorization scales does not differ significantly from the case of the HQET. As an example, for m H = 120 GeV and p H T = 20 GeV, we find for our POWHEG implementation an uncertainty band of [−16%, +21%] around the central value when the factorization and renormalization scales are varied in a range between 0.5 m H and 2 m H while keeping the ratio of the two scales in the range [0. 5,2]. In the case of the HQET the corresponding uncertainty band is [−16%, +20%].
The phenomenological relevance of the finite quark mass effects can be established by comparing the size of the latter to the most accurate estimate of the theoretical uncertainty band for the Higgs transverse momentum distribution. This band can be found e.g. in ref. [43], where NNLO-QCD results in the HQET have been matched analytically with the effects of soft-gluon resummation. For m H = 165 GeV and 10 GeV ≤ p H T ≤ 60 GeV, the band width, for a variation of the factorization, renormalization and resummation scales in a range between 0.5 m H and 2 m H (while keeping the ratio of any two scales in the range [0.5, 2]), amounts to [−8%, +12%] with respect to the central value computed with all the scales set equal to m H . For the same value of m H and range of p H T , we find that the mass effects on the shape of the Higgs transverse momentum distribution range between −8% and +8%, i.e. they are of the same size as the uncertainty band of ref. [43].

MSSM results
We extended our implementation of the gluon-fusion Higgs production process in POWHEG to describe also the production of the neutral CP-even bosons of the MSSM, h and H. Differently from the case of the SM, the CP-even Higgs boson masses of the MSSM can be predicted in terms of the other parameters of the model. At the tree level it is customary to consider the pseudoscalar mass m A and tan β as free parameters which, together with m Z , determine the masses m h and m H of the CP-even Higgs bosons and their couplings to the quarks (as well as to the squarks). Radiative corrections, however, induce in the Higgs masses and couplings a dependence on all of the MSSM parameters (see, e.g., ref. [60]). When studying Higgs boson production in the MSSM, therefore, it is necessary to compute the entire spectrum of masses and couplings of the model in a consistent way, starting from a given set of input parameters. In our numerical analysis, we use the code SoftSusy [61] to compute the MSSM spectrum starting from a set of running parameters expressed in the DR renormalization scheme. However, it is in principle possible to interface our calculation of the Higgs production cross section with other spectrum calculators (such as, e.g., FeynHiggs [62]) that adopt different choices of renormalization scheme for the input parameters.

Modifications in POWHEG
In order to describe the production of the neutral CP-even bosons of the MSSM, three modifications to the SM implementation of POWHEG are required: a rescaling of the Higgs-quark couplings in the top and bottom quark contributions, the introduction of all the contributions from diagrams involving superpartners, and, finally, a rearrangement of the EW corrections (which have been computed only in the SM case).
For the production of the lightest scalar, h, the normalization factors for the Higgs-quark couplings entering the functions H 1ℓ , H 2ℓ , A gg , A qg , A qq in section 3.1 become where α is the effective mixing angle that diagonalizes the radiatively corrected mass matrix in the CPeven Higgs sector. The corresponding normalization factors for the production of the heaviest scalar, H, can be obtained from eq. (31) through the replacements cos α → sin α in λ t and − sin α → cos α in λ b . Diagrams with a squark running in the loop give an additional contribution to the one-loop form factor H 1ℓ : where the sum runs over the six flavors for the squarksq i and the two mass eigenstates (i = 1, 2) for each flavor. The variables yq i and xq i are defined in analogy to y q and x q in eq. (23), and the couplings for the stop and sbottom mass eigenstatest 1 andb 1 to the lightest scalar h read In the equations above µ is the higgsino mass parameter in the MSSM superpotential, A q (for q = t, b) are the soft SUSY-breaking Higgs-squark couplings, θ q are the left-right squark mixing angles and θ W is the Weinberg angle. The couplings for the squark mass eigenstatest 2 andb 2 can be obtained from the corresponding couplings fort 1 andb 1 through the replacements sin 2θ q → − sin 2θ q and cos 2θ q → − cos 2θ q . The squark couplings to the heaviest scalar H can be obtained from the squark couplings to h through the replacements − sin α → cos α and cos α → sin α.
The couplings for the up-type and down-type squarks of the first two generations can be obtained from the stop and sbottom couplings, respectively, by setting the quark mass and the squark mixing angle to zero. However, it can be seen from eqs. (32)-(34) that all contributions from the first two generations of squarks are suppressed by the ratio m 2 Z /m 2 q i . Furthermore, there are significant cancellations among the contributions of the four squarks in each generation (indeed, the total contribution vanishes for degenerate squark masses). Therefore, in what follows we neglect the first two generations, and focus on the stop and sbottom contributions.
Additional contributions to the two-loop form factor H 2ℓ arise from diagrams with squarks and gluons, with four squarks, and with quarks, squarks and gluinos. In our POWHEG implementation we use the results of ref. [27] for the stop contributions, obtained in the limit of vanishing Higgs mass, and the results of ref. [31] for the sbottom contributions, obtained via an asymptotic expansion in the superparticle masses.
The functions A 2 , A 4 and A qq entering the real emission contributions in section 3.1 also receive additional contributions from diagrams with a squark running in the loop: where sq i , tq i and uq i are defined in analogy to s q , t q and u q in eq. (28). Explicit expressions for the functions b 0 (sq i , tq i , uq i ), c 0 (sq i , tq i , uq i ) and d 0 (sq i , tq i , uq i ) are given in ref. [10].
Finally, we need to adapt the electroweak correction δ EW to the case of the MSSM. Although a calculation of the contributions to δ EW from diagrams involving superpartners is not currently available, we can obtain a partial estimate of the EW corrections in the MSSM by introducing in the SM result appropriate rescaling factors for the couplings of the Higgs boson. In particular, when the Higgs boson mass is below the threshold for real top production the EW correction in the SM is dominated by the contribution of two-loop diagrams involving light quarks, in which the Higgs boson couples to a gauge boson [17,18,19]. We can therefore approximate the EW correction for the production of the lightest scalar h as where the one-loop form factor H 1ℓ is computed in the MSSM (i.e., it contains both the quark and squark contributions) and the explicit expression for the two-loop EW light-fermion contribution G 2ℓ lf can be found in ref. [19]. In the case of the production of the heaviest scalar H the factor sin(β − α), which rescales the Higgs-gauge boson couplings, must be replaced by cos(β − α). However, we recall that the approximation of including only the light-fermion contributions becomes less justified when m H 2 m t .

MSSM: numerical results
In this section we present numerical results for the production of the lightest CP-even Higgs boson, h, in a representative region of the MSSM parameter space. Events are generated with our implementation of POWHEG, then matched with the PYTHIA PS. We compute the total inclusive cross section, as well as the transverse momentum distribution, for the production of a light Higgs in gluon fusion, and we compare them with the corresponding quantities computed for a SM Higgs boson with the same mass.
where: m Q , m U and m D are the soft SUSY-breaking mass terms for stop and sbottom squarks; X t ≡ A t − µ cot β is the left-right mixing term in the stop mass matrix; M i (for i = 1, 2, 3) are the soft SUSY-breaking gaugino masses. We consider the input parameters in eq. (39) as expressed in the DR renormalization scheme, at a reference scale Q of the order of the squark masses (in particular, we take Q = 500 GeV). Our choice for X t is modeled on the so-called "m max h scenario", in which the stopinduced radiative corrections maximize the mass of the lightest scalar h, allowing it to satisfy the lower bounds from LEP even for relatively low values of the stop masses (for our choices of parameters the physical masses of the two stops are around 280 GeV and 660 GeV, respectively). We consider both signs for the parameter µ, keeping in mind that, in our conventions, the tan β-dependent corrections to the relation between the bottom mass and the bottom Yukawa coupling [63] enhance the Higgs couplings to bottom and sbottoms for µ < 0 and suppress them for µ > 0.
We perform a scan on the parameters that determine the Higgs boson masses and mixing at tree level, m A and tan β, varying them in the ranges 3 90 GeV ≤ m A ≤ 200 GeV and 2 ≤ tan β ≤ 50. For each value of tan β we derive the soft SUSY-breaking Higgs-stop coupling A t from the condition on X t , then we fix the corresponding Higgs-sbottom coupling as For each point in the parameter space, we use the code SoftSusy [61] to compute the physical (i.e., radiatively corrected) Higgs boson masses m h and m H , and the effective Higgs mixing angle α. We obtain from SoftSusy 4 also the MSSM running quark masses m t and m b , expressed in the DR scheme at the scale Q = 500 GeV. The running quark masses are used both in the calculation of the running stop and sbottom masses and mixing angles, and in the calculation of the top and bottom contributions to the form factors for Higgs boson production (the latter are computed using the DR results presented in refs. [27,31]). As discussed in ref. [31], the use of m b (Q) in the one-loop form factor for gluon fusion, H 1ℓ , induces potentially large contributions, enhanced either by tan β or by ln(m 2 b /Q 2 ), in the two-loop form factor H 2ℓ . We checked that our results are not significantly altered if we compute H 1ℓ in terms of the running bottom mass expressed at the lower scale Q = m h .
In figure 6 we plot the ratio of the cross section for the production of the lightest scalar h in the MSSM over the cross section for the production of a SM Higgs boson with the same mass. For a consistent comparison, we adopt the DR scheme in both the MSSM and the SM calculations. The plot on the left is obtained with µ > 0, while the plot on the right is obtained with µ < 0. In order to interpret the plots, it is useful to recall that for small values of m A it is the heaviest scalar H that has SM-like couplings to fermions, while the coupling of h to top (bottom) quarks is suppressed (enhanced) by tan β. In the lower-left region of the plots, with small m A and moderate tan β, the enhancement of the bottom contribution does not compensate for the suppression of the top contribution, and the MSSM cross section is smaller than the corresponding SM cross section. On the other hand, for sufficiently large tan β (in the lower-right region of the plots) the enhancement of the bottom contribution prevails, and the MSSM cross section becomes larger than the corresponding SM cross section. For µ < 0 the coupling of h to bottom quarks is further enhanced by the tan β-dependent threshold corrections [63], and the ratio between the MSSM and SM predictions can significantly exceed a factor of ten. Finally, for sufficiently large m A , i.e. when the couplings of h to quarks approach their SM values, the MSSM cross section is smaller than the SM cross section.
It is interesting to note that for intermediate values of m A there is a band along which the two cross sections are similar to each other. Indeed, the observation of a scalar particle with cross section in agreement with the SM prediction does not necessarily imply that the Higgs boson is the SM one. However, as will be discussed below, a more detailed study of the Higgs kinematic distributions can help discriminate between the two models.
To assess the genuine effect of the squark contributions (as opposed to the effect of the modifications  in the Higgs-quark couplings), we plot in figure 7 the ratio of the full MSSM cross section for h production over the approximated MSSM cross section computed with only quarks running in the loops. As in figure 6, the plot on the left is obtained with µ > 0, while the plot on the right is obtained with µ < 0. We observe that, in most of the considered region of the MSSM parameter space, the squark contributions reduce the total cross section. We identify three regions: i) for sufficiently large tan β and sufficiently small m A the squark contribution is modest, ranging between −10% and +5%; this region roughly coincides with the one in which the total MSSM cross section is dominated by the tan β-enhanced bottom quark contribution, and is larger than the SM cross section; ii) a transition region, where the corrections rapidly become as large as −30%; this region coincides with the one in which the SM and MSSM cross sections are similar to each other; iii) for sufficiently large m A the squark correction is almost constant, ranging between −40% and −30%; this region coincides with the one in which the MSSM cross section is smaller than the corresponding SM cross section. We now discuss the distribution of the transverse momentum p h T of a light scalar h, considering two distinct scenarios. First, we take a point in the MSSM parameter space (m A = 200 GeV, tan β = 10 and µ > 0) in which the coupling of h to the bottom quark is not particularly enhanced with respect to the SM value, so that the bottom contribution to the cross section is not particularly relevant. Because a light Higgs boson cannot resolve the top and squark vertices, unless we consider very large transverse momentum, we expect the form of the p h T distribution to be very similar to the one for a SM Higgs boson of equal mass, the two distributions just differing by a scaling factor related to the total cross section. This is illustrated in the left plot of figure 8, where we show the ratio of the transverse momentum distribution for h over the transverse momentum distribution for a SM Higgs boson of equal mass. In the right plot of figure 8 we show the ratio of the corresponding shapes, i.e. the distributions normalized to the corresponding cross sections. This ratio, as expected, is close to one in most of the p h T range.    We then consider the opposite situation, namely when the coupling of h to the bottom quark is significantly enhanced. In this situation two tree-level channels, i.e. bb → gh and bg → bh, can also contribute to the production mechanism and influence the shape of the p h T distribution [37]. Leaving a study of the effects of those additional channels to a future analysis, we will now illustrate how the kinematic distribution of the Higgs boson can help discriminate between the SM and the MSSM. The six plots in figure 9 correspond to different points in the (m A , tan β) plane characterized by the fact that the MSSM and SM predictions for the total cross section agree with each other within 5% (therefore, we are effectively comparing the shape of the transverse momentum distributions). The three plots on the left are obtained with µ > 0, while the three plots on the right are obtained with µ < 0. The Higgs boson masses corresponding to these points range between 114 and 122 GeV (i.e., a SM Higgs with the same mass as h would not yet be excluded by direct searches). Figure 9 shows that the region at small p h T receives a positive correction with respect to the SM result for moderate values of tan β. The correction decreases with increasing tan β and eventually becomes negative at large tan β for µ < 0. The region at large p h T shows an opposite behavior with respect to tan β. In figure 10 we show, for the same six points in the (m A , tan β) plane as in figure 9, the ratio of the p h T distribution over the approximate distribution computed with only quarks running in the loops. Even though the light Higgs boson cannot resolve the squark loops, we see that the squark contributions can modify the form of the p h T distribution, because of the interference with the bottom contribution. In particular, we observe that the squark contributions may yield a positive correction on the distribution at small p h T , which turns negative for larger values of the transverse momentum. The negative correction becomes quite flat for p h T > 100 GeV, and in the µ < 0 case it reaches a −50% effect at very large p h T .

Conclusions
We have presented a new implementation 5 in the POWHEG approach of the process of Higgs boson production via gluon fusion in the SM and in the MSSM. In the NLO-QCD contributions, we have retained the exact dependence on all the particle masses in the one-loop diagrams with real-parton emission and in the two-loop diagrams with quarks and gluons, whereas we have employed the approximation of vanishing Higgs mass in the two-loop diagrams involving superpartners. We have also included the effects due to the two-loop EW corrections. The exact mass dependence at NLO QCD and its matching with the multiple gluon emission produced by the PYTHIA PS have important effects on the total and differential cross sections of the Higgs boson. In the SM, the exact dependence on the bottom-quark mass induces, for a light Higgs boson, a non-trivial distortion in the shape of the transverse momentum distribution in the small-p H T region. The effect is comparable in size with the current estimate of the theoretical uncertainty on this observable, which was derived in the HQET limit in ref. [43]. In the case of a heavy Higgs boson, the role of the bottom quark is negligible, but the exact dependence on the top-quark mass yields, with respect to the HQET results, an increase of the transverse momentum distribution at small p H T , and a large negative correction at large p H T . In the MSSM, our code allows for a systematic study of the parameter space of the model in a realistic experimental setup. As an illustration, we considered representative choices in the MSSM parameter space, modeled on the so-called m max h scenario. We studied the role of the bottom diagrams and the impact of the inclusion of diagrams involving superpartners at NLO QCD, both on the total and on the differential cross sections. In the large-tan β regime, where the role of the bottom quark is very relevant, the differential distributions can receive large corrections, which cannot be described in the HQET approximation. A detailed study of the Higgs kinematic distributions could help discriminate between the SM and the MSSM, in case a scalar particle with a cross section compatible with the SM prediction is observed.