Towards precise predictions for Higgs-boson production in the MSSM

We study the production of scalar and pseudoscalar Higgs bosons via gluon fusion and bottom-quark annihilation in the MSSM. Relying on the NNLO-QCD calculation implemented in the public code SusHi, we provide precise predictions for the Higgs-production cross section in six benchmark scenarios compatible with the LHC searches. We also provide a detailed discussion of the sources of theoretical uncertainty in our calculation. We examine the dependence of the cross section on the renormalization and factorization scales, on the precise definition of the Higgs-bottom coupling and on the choice of PDFs, as well as the uncertainties associated to our incomplete knowledge of the SUSY contributions through NNLO. In particular, a potentially large uncertainty originates from uncomputed higher-order QCD corrections to the bottom-quark contributions to gluon fusion.


Introduction
The recent discovery of a Higgs boson with mass around 125.5 GeV by the ATLAS and CMS experiments at the Large Hadron Collider (LHC) [1,2] puts new emphasis on the need for precise theoretical predictions for Higgs production and decay rates, both in the Standard Model (SM) and in plausible extensions of the latter such as the Minimal Supersymmetric Standard Model (MSSM). The current status of these calculations is summarized in the reports of the LHC Higgs Cross Section Working Group (LHC-HXSWG) [3,4,5].
In the SM, the main mechanism for Higgs production at hadron colliders is gluon fusion [6], where the coupling of the gluons to the Higgs is mediated by loops of heavy quarks, primarily top and bottom. The knowledge of this process includes: the next-to-leading order (NLO) QCD contributions [7] computed for arbitrary values of the Higgs and quark masses [8,9,10,11]; the next-to-next-toleading order (NNLO) QCD contributions due to top-quark loops, in the heavy-top limit [12,13] and including finite top-mass effects [14]; soft-gluon resummation effects [15] and estimates of the nextto-next-to-next-to-leading order (NNNLO) QCD contributions [16]; the first-order electroweak (EW) contributions [17,18,19,20,21] and estimates of the mixed QCD-EW contributions [22].
The Higgs sector of the 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 scalars, h and H, one neutral pseudoscalar, A, and two charged scalars, 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, one of the most important production mechanisms for the neutral Higgs bosons is gluon fusion, mediated by loops involving the top and bottom quarks and their superpartners, the stop and sbottom squarks. However, for intermediate to large values of tan β bottom-quark annihilation can become the dominant production mechanism for the neutral Higgs bosons that have enhanced couplings to down-type fermions.
If the third-generation squarks have masses around one TeV or even larger, their contributions to the gluon-fusion process are suppressed, and a sufficiently accurate determination of the cross section can be achieved by rescaling the SM results for the top-and bottom-quark contributions by appropriate Higgs-quark effective couplings. If, on the other hand, some of the squarks have masses of the order of a few hundred GeV -a scenario not yet excluded by the direct searches at the LHC -a precise calculation of the contributions to the gluon-fusion cross section from diagrams involving squarks becomes mandatory. The NLO-QCD contributions to scalar production arising from diagrams with colored scalars and gluons were first computed in the vanishing-Higgs-mass limit (VHML) in ref. [23], and the full Higgs-mass dependence was included in later calculations [10,11,24]. For what concerns pseudoscalar production, the NLO-QCD contributions arising from diagrams with quarks and gluons are known [8,9,10,11] while diagrams involving only squarks and gluons do not contribute to the gluon-fusion process due to the structure of the pseudoscalar couplings to squarks. In contrast, a full calculation of the contributions to either scalar or pseudoscalar production arising from two-loop diagrams with quarks, squarks and gluinos -which can involve up to five different particle masses -is still missing. Calculations based on a combination of analytic and numerical methods were presented in refs. [25,26], but neither explicit analytic formulae nor public computer codes implementing the results of those calculations have been made available so far.
Approximate results for the quark-squark-gluino contributions can however be obtained assuming the presence of some hierarchy between the Higgs mass and the masses of the particles running in the loops. If the Higgs boson is lighter than all the particles in the loops, it is possible to expand the result in powers of the Higgs mass, with the first term in the expansion corresponding to the VHML. This limit was adopted in refs. [27,28,29] for the calculation of the top-stop-gluino contributions to scalar production and in refs. [30,31] for the analogous calculation of pseudoscalar production. Refs. [29,31] also discussed the reliability of the VHML by considering the next term in the expansion in the Higgs mass.
While an expansion in the Higgs mass is a viable approximation in the computation of the topstop-gluino contributions to the production of the lightest scalar h, it might not be applicable to the production of the heaviest scalar H and of the pseudoscalar A, if their mass is comparable to the mass of the top quark. Moreover, an expansion in the Higgs mass is certainly useless in the calculation of the bottom-sbottom-gluino contributions, due to the presence of a light bottom quark. All of these limitations can, however, be overcome with an expansion in inverse powers of the superparticle masses. Since it does not assume any hierarchy between the Higgs mass and the mass of the quark in the loop, such an expansion is applicable to both top-stop-gluino and bottom-sbottom-gluino contributions, as long as the squarks and the gluino are heavier than the considered Higgs boson and the top quark. Results for scalar production based on an expansion in the superparticle masses were presented in refs. [32,33,34], and analogous results for pseudoscalar production were presented in ref. [31].
In order to improve the accuracy of the MSSM prediction for the gluon-fusion cross section, and to allow for a meaningful comparison with the SM prediction, several contributions beyond the NLO in QCD should be included. The NNLO-QCD contributions to scalar production arising from diagrams with top quarks and the subset of EW contributions arising from diagrams with light quarks can be obtained from the corresponding SM results with an appropriate rescaling of the Higgs couplings to quarks and to gauge bosons. The NNLO-QCD top-quark contributions to pseudoscalar production have also been computed [35]. Approximate results beyond the NLO in QCD also exist for the contributions of diagrams involving superparticles. A first estimate of the NNLO-QCD contributions of diagrams involving stop squarks was presented in ref. [36], and an approximate calculation of those contributions, assuming the VHML and specific hierarchies among the superparticle masses, was recently presented in refs. [37,38]. Furthermore, a subset of potentially large tan β-enhanced contributions from diagrams involving sbottom-gluino or stop-chargino loops can be resummed in the LO cross section by means of an effective Higgs-bottom coupling [39,40,41].
In a significant part of the MSSM parameter space, the couplings of the heavier neutral Higgs bosons H and A to bottom quarks are enhanced by tan β with respect to the corresponding coupling of the SM Higgs, while their couplings to top quarks are suppressed by tan β. When that is the case, the bottom-quark contributions to the gluon-fusion process -which for a SM-like Higgs with mass around 125.5 GeV amount to roughly 7% of the cross section -can dominate over the top-quark contributions. The bottom-quark contributions are subject to large QCD corrections enhanced by powers of ln(m 2 φ /m 2 b ), where φ denotes a generic Higgs boson, and so far they have been computed only at the NLO [8,9,10,11]. As a result, the uncomputed higher-order QCD corrections to the bottom-quark contributions can become the dominant source of uncertainty in the cross section for the production of heavy MSSM Higgs bosons in gluon fusion.
On the other hand, as mentioned earlier, when the couplings to bottom quarks are sufficiently enhanced the production of MSSM Higgs bosons through bottom-quark annihilation dominates over gluon fusion. In the four-flavor scheme (4FS), where one does not consider the bottom quarks as partons in the proton, the process is initiated by two gluons or by a light quark-antiquark pair, and the cross section is known at the NLO in QCD [42]. In the five-flavor scheme (5FS), where the bottom quarks are in the initial partonic state, the cross section is known up to the NNLO in QCD [43,44]. The use of bottom-quark parton density functions (PDFs) in the 5FS allows to resum terms enhanced by ln(m 2 φ /m 2 b ) that would arise in the 4FS when one or both bottom quarks are collinear to the incoming partons. As in the case of gluon fusion, the tan β-enhanced contributions from diagrams involving superpartners can be resummed in the LO result by means of an effective Higgs-bottom coupling. The remaining one-loop contributions from superpartners have been found to be small [45].
A considerable effort has been devoted over the years to making the existing calculations of Higgs production available to the physics community in the form of public computer codes. In the case of the SM, NNLO-QCD predictions of the total cross section for gluon fusion, including various refinements such as EW corrections and finite top-mass effects, are provided, e.g., by HIGLU [46], ggh@nnlo [47], HNNLO [48] and iHixs [49]. The code bbh@nnlo [50] provides instead a NNLO-QCD prediction of the total cross section for Higgs production in bottom-quark annihilation in the 5FS. For what concerns the production of MSSM Higgs bosons via gluon fusion, HIGLU implements the results of ref. [24] for the NLO-QCD contributions arising from diagrams with squarks and gluons, as well as the results of refs. [40,41] for the resummation of the tan β-enhanced squark contributions in an effective Higgsbottom coupling.
More recently, two codes that compute the cross section for Higgs production including approximate results for the contributions of diagrams with quarks, squarks and gluinos have become available. As described in ref. [51], the NLO-QCD [11,29,31,32,34] and EW [18,21] contributions to Higgsboson production via gluon fusion in the SM and in the MSSM have been implemented in a module for the so-called POWHEG BOX [52], a framework for consistently matching NLO-QCD computations of matrix elements with parton-shower Monte Carlo generators, avoiding double counting and preserving the NLO accuracy of the calculation. The code SusHi [53], on the other hand, computes the cross section for Higgs-boson production in both gluon fusion and bottom-quark annihilation, in the SM and in the MSSM. In the case of gluon fusion, SusHi includes the exact results of ref. [9] for the NLO-QCD contributions of two-loop diagrams with top and bottom quarks, and the approximate results of refs. [28,31,34] and refs. [31,32] for the NLO-QCD contributions of two-loop diagrams with stop and sbottom squarks, respectively. The NLO-QCD contributions of one-loop diagrams with emission of an additional parton are taken from ref. [33]. The NNLO-QCD contributions from diagrams with top quarks are included via a call to ggh@nnlo, and the corresponding contributions from diagrams with stop squarks are estimated following ref. [36]. Finally, the known SM results for the EW contributions [18,20,21] are adapted to the MSSM by rescaling the Higgs couplings to top quarks and to gauge bosons. In the case of bottom-quark annihilation, SusHi obtains from bbh@nnlo the NNLO-QCD result valid in the SM, then rescales it by an effective Higgs-bottom coupling that accounts for the tan β-enhanced squark contributions [39,40].
In this paper we use SusHi for a precise study of scalar and pseudoscalar Higgs production in the MSSM. In section 2 we present predictions for the total inclusive cross section for Higgs production in six benchmark scenarios compatible with the LHC results, focusing in particular on a scenario with relatively light stops where the effect of the SUSY contributions can be significant. In section 3 we provide a detailed discussion of the sources of theoretical uncertainty in the calculation of the total cross section for Higgs-boson production in the MSSM. We examine the dependence of the cross sections for gluon fusion and bottom-quark annihilation on the renormalization and factorization scales, on the precise definition of the Higgs-bottom coupling and on the choice of PDFs, as well as the uncertainty associated to our incomplete knowledge of the SUSY contributions through NNLO. In particular, we point out a potentially large uncertainty arising from uncomputed higher-order QCD corrections to the bottom-quark contributions to gluon fusion, which can affect the interpretation of the searches for the MSSM Higgs bosons in scenarios where their couplings to bottom quarks are enhanced with respect to the SM. In section 4 we present our conclusions. Finally, in the appendix we list the cross sections and uncertainties for the production of the three neutral Higgs bosons in selected points of the parameter space for the six benchmark scenarios.

Higgs-boson production in viable MSSM scenarios
The discovery of a neutral scalar with mass around 125.5 GeV puts the studies of the Higgs sector of the MSSM in an entirely new perspective. In order to remain viable, a point in the MSSM parameter space must now not only pass all the experimental bounds on superparticle masses, but also lead to the prediction of a scalar with mass, production cross section and decay rates compatible with those measured at the LHC. In particular, the relatively large mass of the SM-like scalar discovered at the LHC implies either stop masses of the order of 3 TeV -which would result in a negligible stop contribution to the production cross section -or a large value of the left-right mixing term in the stop mass matrix (see, e.g., refs. [54,55]). In the latter case, at least one of the stops could have a mass as low as a few hundred GeV, and induce a significant contribution to the gluon-fusion cross section.  Table 1: Choices of MSSM parameters for the benchmark scenarios proposed in ref. [56].
FeynHiggs takes by default equal to m t , while m A is identified with the pole mass of the pseudoscalar. Finally, the choice of renormalization scheme for mg amounts to a higher-order effect, because the gluino mass enters only the two-loop part of the corrections. A detailed description of the six benchmark scenarios adopted in our analysis can be found in the paper where they were originally proposed, ref. [56]. All of the scenarios are characterized by relatively large values of the ratio X t /M S , ensuring that the mass of the SM-like Higgs falls within the required range without the need for extremely heavy stops. In addition, the masses of the gluino and of the first-two-generation squarks are set to 1.5 TeV, large enough to evade the current ATLAS [71,72] and CMS [73,74,75] bounds. The prescriptions of ref. [56] for the parameters M S , X t , µ and for the soft SUSY-breaking wino mass M 2 are listed in table 1. We vary the parameters tan β and m A within the ranges 2 ≤ tan β ≤ 50 , In all scenarios the Higgs-sbottom-sbottom coupling A b is set equal to A t , the left-right mixing of the first-two-generation squarks is neglected and the bino mass M 1 is obtained from the GUT relation , with the exception of the fourth scenario where we set M 1 = 340 GeV. 3 Finally, the choices of ref. [56] for the soft SUSY-breaking parameters in the slepton sector have a very small impact on the predictions for the Higgs masses and production cross sections, therefore we do not report them here.
The fourth scenario in table 1, denoted as light stop, deserves a special discussion. In this scenario the two stop masses are 324 GeV and 672 GeV; the sbottom masses depend on tan β, but the lightest sbottom is always heavier than 450 GeV, while the heaviest one is always lighter than 550 GeV. With such relatively low masses, loops involving squarks can give a sizable contribution to the cross section for Higgs production, but we have to worry about the exclusion bounds from the LHC. Indeed, the ATLAS and CMS collaborations have presented preliminary results for the searches of direct stopand sbottom-pair production, based on the full 8-TeV data sample, considering the decay chains The allowed values of the stop and sbottom masses depend on the chargino and neutralino masses, as well as on the branching ratios for the different squark decays. With the choice of parameters in table 1, M 2 = µ = 400 GeV, together with M 1 = 340 GeV, the masses for the lightest chargino and neutralino have a mild dependence on tan β, but they stay within the ranges m χ ± 1 ≈ 341 -346 GeV and m χ 0 1 ≈ 316 -320 GeV for tan β > 10. In this case the lightest stop decays almost entirely through the loop-induced, flavor-violating channelt 1 → c χ 0 1 . This channel has been investigated by ATLAS [78] and CMS [79], but the resulting bounds only reach to values of mt 1 around 250 GeV. For the lightest sbottom, the two-body decaysb 1 →t 1 W andb 1 → b χ 0 j (with j up to 3 or 4) are kinematically open. The direct decay ofb 1 to the lightest neutralino would be constrained by the searches in refs. [73,80], but i) that channel is never dominant in the considered range of parameters and ii) the experimental bounds only reach to values of m χ 0 1 below 280 GeV. Finally, the heaviest stop and sbottom can decay through a multitude of channels, and their direct decays to χ 0 1 or χ ± 1 are significantly suppressed.

Cross section for Higgs production
We are now ready to present our precise predictions for the production of MSSM Higgs bosons at the LHC. As mentioned earlier, we rely on the code SusHi, 4 which includes all of the available NLO-QCD contributions to the gluon-fusion process, supplemented with the known SM results for the NNLO-QCD contributions in the heavy-top limit and for the EW contributions (both adapted to the MSSM by appropriately rescaling the Higgs couplings). While the results implemented in SusHi for the NNLO-QCD top contributions are strictly valid only for a Higgs mass below the top threshold, m φ < 2 m t , a comparison with the NLO results suggests that they provide a decent approximation also for larger values of the Higgs mass [81,82]. The NNLO-QCD contributions from stop loops are estimated following ref. [36], i.e., neglecting the contributions of three-loop diagrams but retaining the NNLO contributions that arise from the product of lower-order terms. We have also checked that, when all of the NNLO-QCD contributions are omitted, the results of SusHi for the gluon-fusion cross section agree with those of the calculation implemented in the POWHEG BOX [51], which includes the same NLO-QCD and EW contributions. For what concerns the bottom-quark annihilation process, SusHi includes the NNLO-QCD results valid in the SM within the 5FS, also rescaled by the effective Higgs-bottom couplings of the MSSM.
In our study, we fix the center-of-mass energy of the proton-proton collisions to 8 TeV. While the numerical value of the total cross section for Higgs production does obviously depend on the collision energy, we have checked that the relative importance of the various contributions to the production processes and their qualitative behavior over the MSSM parameter space do not change substantially if we set the energy to 13 TeV. By default, we use the MSTW2008 set of PDFs [83], and we fix the renormalization and factorization scales entering the gluon-fusion cross section to µ R = µ F = m φ /2 [13,84], where φ = {h, H, A} denotes the considered Higgs boson. For bottomquark annihilation, the central values of the scales are chosen as µ R = m φ and µ F = m φ /4 [43,44,85]. In the calculation of the gluon-fusion cross section we relate the bottom Yukawa coupling to the pole mass M b , computed at the three-loop level [86] from the input value for the running mass, m b (m b ). In the case of bottom-quark annihilation, on the other hand, we relate the bottom Yukawa coupling to m b (m φ ), in turn obtained from m b (m b ) via four-loop renormalization-group evolution [87]. In both cases, the tan β-enhanced SUSY corrections to the relation between mass and Yukawa coupling of the bottom quark are included following refs. [39,40]. The theoretical uncertainties associated to the choice of PDFs, to the variation of the renormalization and factorization scales and to the definition of the bottom Yukawa coupling will be discussed in detail in section 3.
In figures 1 and 2 we show the total cross section -i.e., the sum of gluon fusion and bottom-quark annihilation -for the production of the scalars (h, H) and of the pseudoscalar (A), respectively, as contour plots in the m A -tan β plane. For the other MSSM parameters, we adopt the light-stop scenario described in section 2.1. Tables for the numerical values of the cross section (and the corresponding uncertainties) in all of the six benchmark scenarios are given in the appendix. In the two plots of figure 1, referring to h (left) and H (right) production, the red lines are contours of equal mass for the corresponding scalar. In this scenario, the prediction for the mass of the lightest scalar reaches a maximum of 123.8 GeV at large tan β. The heaviest-scalar mass grows with m A , and we show only the contour corresponding to 126 GeV to avoid clutter (for large m A , the contours are roughly at m H ≈ m A and independent of tan β). The x-axis of the plot for h production ends at m A = 300 GeV because, for larger values, the cross section becomes essentially independent of m A . The x-axis of the plots for H and A ends at m A = 500 GeV because the expansion in the SUSY masses used to approximate the two-loop squark contributions in SusHi becomes unreliable when the Higgs mass approaches the lowest squark-mass threshold, which in the light-stop scenario corresponds to 2 mt 1 ≈ 650 GeV. The theoretical uncertainty associated with this approximation will be discussed in section 3.4.
The qualitative behavior of the cross sections in figures 1 and 2 can be easily interpreted considering the relations between the scalar and pseudoscalar masses in the MSSM Higgs sector, and how each of the Higgs bosons couples to the top and bottom quarks (the squark contributions are generally sub-dominant, as will be discussed below). In the so-called decoupling limit, m A m Z , the lightest scalar h has SM-like couplings to quarks, while its mass is essentially independent of m A and, for tan β > ∼ 10, depends only weakly on tan β. The cross section for h production (left plot in figure 1) varies very little in this region, and differs from the SM result for a Higgs boson of equal mass only because of the squark contributions to the gluon-fusion process. For m A < ∼ 130 GeV, on the other hand, the couplings of h to top (bottom) quarks are non-standard, being suppressed (enhanced) by tan β. In this narrow region the total cross section for h production is dominated by the contributions of the diagrams that involve the Higgs-bottom coupling, and it grows significantly with tan β.  The behavior of the cross section for H production in the m A -tan β plane (right plot in figure 1) is different from -and somewhat complementary to -the one for h production. In the strip where m A < ∼ 130 GeV, the heaviest scalar has a mass around 125 GeV and significant couplings to both top and bottom quarks, and the cross section for its production grows with tan β. For larger m A , on the other hand, m H grows together with m A , and the couplings of H to top (bottom) quarks are suppressed (enhanced) by tan β. The total cross section for H production is therefore dominated, already for moderate tan β, by the contributions of the diagrams that involve the Higgs-bottom coupling. The latter grow significantly with tan β, but decrease with m A , being suppressed by powers of the ratio m 2 b /m 2 H . Finally, the pseudoscalar couplings to top (bottom) quarks are suppressed (enhanced) by tan β for all values of m A . Therefore, the behavior of the cross section for A production in the m Atan β plane, see figure 2, resembles the behavior of h production when m A < ∼ 130 GeV, and the one of H production for larger m A : in both cases, the cross section grows with tan β, but decreases with m A .
To disentangle the effects of the two main production channels for the MSSM Higgs bosons, we show in figures 3 and 4 the ratio between the gluon-fusion cross section and the sum of gluon-fusion and bottom-quark-annihilation cross sections in the light-stop scenario, again as contour plots in the m Atan β plane. Predictably, the plots reflect the behavior of the coupling of the considered Higgs boson to bottom quarks. The left plot in figure 3 shows that, when m A is large enough that the couplings of the lightest scalar are SM-like, gluon fusion is by far the dominant process for h production, and the contribution of bottom-quark annihilation amounts only to a few percent. Only in the strip with m A < ∼ 130 GeV and tan β > ∼ 8, where the coupling of h to bottom quarks is sufficiently enhanced by tan β, does bottom-quark annihilation become the dominant process. Conversely, bottom-quark annihilation gives the largest contribution to the cross section for H production (right plot in figure 3) when m A > ∼ 130 GeV and tan β > ∼ 6, while in the case of A production (figure 4) the cross section is dominated by bottom-quark annihilation already for m A > ∼ 100 GeV, as long as tan β > ∼ 5 -8. To assess the relevance of the squark contributions to the gluon-fusion cross section in the light-stop scenario, we show in figures 5 and 6 the ratio of the total gluon-fusion cross section over the cross section computed including only the contributions of quarks (with appropriate rescaling of the Higgsquark couplings). The left plot of figure 5 shows that -in this scenario characterized by relatively light squarks -the interference between the top and stop contributions can reduce the cross section for h production by as much as 20% in the decoupling region with large m A and tan β. Remarkably, in this region the partial NNLO-QCD contributions from stop loops that we include following ref. [36] account by themselves for a 6% suppression of the cross section. The theoretical uncertainty associated to these contributions will be discussed in section 3.4. For what concerns H production (right plot of figure 5), the squark contributions reduce the cross section by up to 30% for low values of m A , and the suppression becomes even stronger with increasing pseudoscalar mass. In particular, near the lower-right corner of the plot, where m A > ∼ 420 GeV and tan β ranges between 6 and 20, the interference between the quark and squark contributions induce a suppression of the cross section by 70 -80%. In this region the top contribution is suppressed by tan β, while the bottom contribution is     suppressed by m 2 b /m 2 H and only moderately enhanced by tan β, so they both become comparable in size with the stop contribution. The resulting gluon-fusion cross section is rather small, of the order of a few femtobarns. Finally, figure 6 shows that, in the case of A production, the effect of the squark contributions on the cross section for gluon fusion in the light-stop scenario is always less than 10%. This is due to the fact that the pseudoscalar couples only to two different squark-mass eigenstates, while gluons couple only to pairs of the same squarks. Therefore, there is no squark contribution to the gluon-fusion process at the LO, and the whole effect in figure 6 arises from two-loop diagrams.
For a SM Higgs boson sufficiently lighter than the top threshold, the EW corrections to gluon fusion are well approximated [20,21] by the contributions of two-loop diagrams in which the Higgs couples to EW gauge bosons, which in turn couple to the gluons via a loop of light quarks (including the bottom). In SusHi, these contributions are incorporated in the MSSM calculation of the gluon-fusion cross section by rescaling the two-loop EW amplitude given in ref. [21] with the appropriate Higgsgauge boson couplings. 5 In figure 7 we investigate the impact of the light-quark EW contributions on the production of the scalars h and H, plotting the ratio of the gluon-fusion cross sections computed with and without those contributions, in the m A -tan β plane for the light-stop scenario. The figure shows that the EW corrections tend to increase the cross section, and their impact depends mainly on the strength of the coupling of the considered scalar to gauge bosons. In the case of h production (left plot) the EW corrections become fairly constant, around 6%, in the region of sufficiently large m A where the lightest scalar has SM-like couplings. Conversely, in the case of H production (right plot) the EW corrections reach a comparable value only in the strip of very low m A , and they quickly drop below 1% as soon as m A > ∼ 150 GeV. On the other hand, since the pseudoscalar does not couple to two gauge bosons at tree level, there are no EW contributions from light-quark loops to its production.
For what concerns the remaining sources of EW corrections to gluon fusion, those arising from two-loop diagrams involving top quarks are known to be small for a SM-like Higgs with mass around 125 GeV [20], while in the case of H and A they are suppressed in most of the parameter space by the small (or vanishing) Higgs couplings to top quarks and to gauge bosons. On the other hand, the EW corrections involving the bottom Yukawa coupling, which have not yet been computed because they are negligible for the SM Higgs, could become relevant for the production of H and A. In addition, a full computation of the EW corrections should include the contributions of diagrams involving superparticles. The non-decoupling SUSY effects that dominate at large tan β are indeed included in an effective Higgs-bottom coupling, as discussed in section 3.2.2, but the remaining contributions, so far uncomputed, could become relevant if some of the superparticles are relatively light.
Results for the Higgs-production cross section in the other benchmark scenarios listed in table 1 can be found in the appendix. In the four scenarios denoted as m max h , m mod+ h , m mod− h and light stau, the couplings of the Higgs bosons to top and bottom quarks and to gauge bosons are rather similar to the ones in the light-stop scenario. Thus, the discussion given above for the qualitative behavior in the m A -tan β plane of the total cross section, of the EW corrections and of the relative importance of gluon fusion and bottom-quark annihilation applies to those four scenarios as well. However, all of the third-generation squarks have masses around 1 TeV, therefore the impact of the SUSY contributions on the gluon-fusion cross section is considerably smaller than in the case of the light-stop scenario. The suppression of the cross section for h production in the decoupling limit never goes beyond 6%. For what concerns H production, the effect of the interference between quark and squark contributions becomes significant only for very large m A and moderate tan β, where the gluon-fusion cross section is tiny anyway. The largest effect, a suppression by 30 -40%, is found in the light-stau scenario for m A > ∼ 850 GeV and 10 < ∼ tan β < ∼ 20, where the cross section is of the order of a tenth of a femtobarn. The SUSY contributions to A production, already small in the light-stop scenario because they only arise at two loops, are further suppressed in the m max h , m mod+ h , m mod− h and light-stau scenarios. In the last scenario in table 1, denoted as tau-phobic, the MSSM parameters are arranged in such a way that, for certain values of m A and tan β, the radiative corrections to the (1, 2) element of the CP-even Higgs mass matrix suppress significantly the mixing angle α, so that the coupling of h to taus -which is proportional to sin α -is in turn suppressed with respect to its SM value. However, the couplings of the scalars to top and bottom quarks are modified as well, in particular the coupling of h to bottom quarks is suppressed. As a result, in the tau-phobic scenario the behavior in the m A -tan β plane of the various contributions to the Higgs-production cross section differs from the one found in the other scenarios. The total cross section for h production shows some enhancement with tan β even for large values of m A , while for small m A the total cross section for H production has a milder dependence on tan β than in the other scenarios. Also, the suppression of the h coupling to bottom quarks makes the contribution of bottom-quark annihilation to h production smaller than in the other scenarios. Finally, the tau-phobic scenario is characterized by third-generation squark masses around 1.5 TeV, and by a value of the superpotential Higgs-mass parameter, µ = 2 TeV, much larger than in the other scenarios. Since µ enters the couplings of the Higgs bosons to squarks, the impact of the SUSY contributions on the cross section for scalar production is -despite the heavier squarks -somewhat larger than in the m max h , m mod+ h , m mod− h and light-stau scenarios, and in the case of pseudoscalar production it is even larger than in the light-stop scenario.

Sources of theoretical uncertainty
Like any other quantity evaluated perturbatively, the cross sections for Higgs production in gluon fusion and bottom-quark annihilation suffer from an intrinsic theoretical uncertainty due to the truncation at finite order in the coupling constants. Typically, the residual dependence on the renormalization and factorization scales is used as an estimate of this uncertainty. In section 3.1 we discuss our study of the scale dependence of the cross sections.
In addition, there are sources of uncertainty that are more specific to the Higgs-production processes considered in this paper. As we discuss in section 3.2, one of the most important sources of uncertainty in the production of Higgs bosons with non-standard couplings to quarks is the dependence of the cross section on the precise definition of the bottom-quark mass and Yukawa coupling. The numerical difference between the pole bottom mass and the running mass computed at a scale of the order of the Higgs mass is more than 40%, and -in a fixed-order calculation of the cross sections -the effect of such a large variation cannot be compensated by the large logarithms that are induced at NLO by counterterm contributions. Furthermore, it is well known that the relation between the bottom mass and the corresponding Yukawa coupling is affected by potentially large, tan β-enhanced SUSY corrections that must be properly resummed. The dependence of the cross sections on the details of the resummation procedure constitutes a further source of uncertainty.
In section 3.3 we discuss the uncertainties associated to the choice of PDF sets. We also investigate the issue of consistency between the pre-defined value of the bottom mass in the PDFs and the value of the mass used to extract the bottom Yukawa coupling.
Finally, in section 3.4 we discuss two sources of uncertainty arising from our incomplete knowledge of the SUSY contributions to gluon fusion. In particular, we assess the validity of the expansion in inverse powers of the SUSY masses used to approximate the contributions of two-loop diagrams involving superparticles. We also estimate the uncertainty associated to the fact that SusHi does not include the contributions of three-loop diagrams involving superparticles.

Scale dependence of the cross section
In this section we study the dependence of the cross section for Higgs production on the renormalization scale µ R at which the relevant couplings in the partonic cross section are expressed, and on the factorization scale µ F entering both the PDFs and the partonic cross section. We recall that, although the complete result for the hadronic cross section does not depend on µ R and µ F , its approximation at a given perturbative order retains a dependence on those scales, which is formally one order higher than the accuracy of the calculation. In a given calculation at fixed order, the two scales are arbitrary, and they are typically fixed at some central valuesμ R andμ F characteristic of the hard scattering process. The variation of the scales around their central values provides an estimate of the size of the uncomputed higher-order contributions.
We discuss separately the cases of gluon fusion (section 3.1.1) and of bottom-quark annihilation (section 3.1.2). In the former, µ R denotes the scale at which we express the strong gauge coupling entering the partonic cross section already at the LO, while in the latter it denotes the scale at which we express both the bottom Yukawa coupling entering at the LO and the strong gauge coupling entering at the NLO. We postpone to section 3.2 a discussion of the dependence of the gluon-fusion cross section on the scale at which we express the bottom Yukawa coupling.

Gluon fusion
The natural hard scale in the production of a Higgs boson φ is obviously of the order of m φ . In our study of gluon fusion we takeμ R =μ F = m φ /2 as central values for the renormalization and factorization scales, because, with this choice, the cross section shows a reduced sensitivity to scale variations and an improved convergence of the perturbative expansion [13]. Moreover, it has been observed that this choice allows to mimic the effects of soft-gluon resummation in the total cross section [84].
We study the impact of the scale variation around the central choice (μ R ,μ F ) following the LHC-HXSWG prescription [3]: we consider seven combinations of renormalization and factorization scales, defined as the set C µ of the pairs (µ R , µ F ) obtainable from the two sets , we treat the variations of the ratio µ R /µ F on the same footing as the variations of the individual scales, discarding the two pairs where the ratio varies by a factor of four around its central value). We then determine the maximal and minimal values of the cross section on the set C µ , and define the relative scale uncertainty of the cross section as    In figures 8 and 9 we show the contours of equal ∆ µ for scalar and pseudoscalar production in the m A -tan β plane, fixing the MSSM parameters as in the light-stop scenario. The qualitative features of the plots can be understood by considering that the top, bottom, SUSY and EW contributions to the gluon-fusion cross section are known at different orders in the perturbative expansion. In particular, the top contribution is included in SusHi with full mass dependence through O(α 3 s ) (i.e., NLO) and in the VHML at O(α 4 s ) (i.e, NNLO). Its residual scale dependence amounts to an O(α 5 s ) effect, with the exception of some mass-dependent effects at O(α 4 s ), which are known to be numerically small [14]. The bottom and sbottom contributions are included at the NLO and they account for an O(α 4 s ) effect. The stop contributions are included through the NNLO, see section 3.4, but their effect on scale dependence is also of O(α 4 s ) because we neglect the genuine three-loop terms. Finally, while the EW corrections are computed at O(αα 2 s ), their inclusion as a fully factorized term at the NLO causes their effect on scale variation to be of O(αα 4 s ), numerically very small. As a consequence of the varying accuracy of the different contributions, the scale uncertainty for the production of a given Higgs boson depends on which contribution plays the dominant role in the considered region of the m A -tan β plane. The uncertainty is lowest, around 10 -20%, where the top contribution dominates: this is the case for h production (left plot in figure 8) in the decoupling region, where the uncertainty stabilizes to roughly 16% at large m A (i.e., slightly smaller than the 18% we obtain for the same Higgs mass in the SM); for H production (right plot in figure 8) in the strip with m A < ∼ 120 GeV, as well as when tan β < ∼ 10 and m A < ∼ 400 GeV; for A production (figure 9) in the strip with tan β < ∼ 10. In contrast, the scale uncertainty exceeds 20% in the regions where the bottom contribution is enhanced or downright dominant: at large tan β for H and A production, and at small m A for h production.
The plots for H and A production in figures 8 and 9 show additional structures. In the case of H production, the scale uncertainty becomes very large for 8 < ∼ tan β < ∼ 16 and m A > ∼ 460 GeV. As appears from figure 5, this region is characterized by a significant cancellation between the top, bottom and stop contributions to the gluon-fusion amplitude, resulting in a very small NLO cross section and an enhanced sensitivity to higher-order effects. In the case of A production, the structure visible for m A ≈ 350 GeV is associated to the cusp-like behavior of the top contribution to the gluonfusion amplitude around the threshold m A = 2 m t . Another feature of H and A production, partially overshadowed by the structures described above, is a tendency towards smaller scale uncertainties for larger pseudoscalar (and hence scalar) masses. This is due to the fact that the strong gauge coupling -which controls the size of the higher-order effects that we are estimating -is evaluated at a scale proportional to the mass of the considered Higgs boson, and gets smaller when the scale increases.
The other scenarios were studied following the same procedure, and the results are qualitatively similar. For h production, the scale dependence in the decoupling region is similar to, or even bigger than, the one in the SM. For H production, due to the different interplay of quark and squark contributions, the cancellations that in the light-stop scenario cause the region of very large uncertainty for 8 < ∼ tan β < ∼ 16 and m A > ∼ 460 GeV occur at higher values of m A . Finally, a study of independent variations of the renormalization and factorization scales shows that, in a large fraction of the parameter space, the former yield a much larger uncertainty than the latter. The factorization-scale uncertainty is smaller in size than the renormalization-scale uncertainty already at the LO, and it is further reduced by the inclusion of higher-order terms.

Bottom-quark annihilation
In SusHi, the cross section for Higgs production in bottom-quark annihilation is implemented at NNLO-QCD in the 5FS. Our default choice for the central scales isμ R = m φ andμ F = m φ /4, following the observation that radiative corrections are particularly small for this value of the factorization scale [43,44,85]. To study the uncertainty associated to the variation of the scales, we consider seven combinations corresponding to all possible pairings of , with the additional constraint that 2 ≤ µ R /µ F ≤ 8 (again, we discard the two pairs with the largest variation of µ R /µ F around its central value, which in this case is 4 ). We then determine the scale uncertainty ∆ µ in analogy to eqs. (2) and (3).
Differently from the case of gluon fusion, the scale uncertainty of bottom-quark annihilation depends very weakly on tan β. This is due to the fact that, in eq. (3), the tan β-dependence of the cross section via the effective Higgs-bottom coupling cancels out in the ratio, leaving only a mild, indirect dependence -only for scalar production -via the value of the Higgs mass that determines µ R and µ F .
In figures 10 and 11 we show the scale dependence of the cross section for scalar and pseudoscalar production, respectively, as a function of m A in the light-stop scenario with tan β = 20. In the upper part of each plot, the solid line denotes the cross section for bottom-quark annihilation computed with the central scale choice (μ R ,μ F ), while the yellow band around the solid line is delimited by the maximal and minimal cross sections σ + and σ − , defined in analogy to eq. (2). The lower part of each plot shows the relative variation of the cross section with respect to the central value (i.e., the total width of the yellow band corresponds to ∆ µ ). While the values of the total cross section do of course depend on the chosen benchmark scenario, the relative scale variation is essentially the same in all scenarios, due to the above-mentioned cancellation of the dependence on the effective Higgs-bottom coupling.
The left plot in figure 10 shows that the relative scale uncertainty of the cross section for h production can be as large as 30% for low values of m A , then it stabilizes to roughly 18% in the decoupling region where m h becomes independent of m A . In contrast, the relative scale uncertainty of the cross section for the production of H (right plot in figure 10) and A (figure 11) decreases as m A (and hence m H ) increases. As already mentioned for the case of gluon fusion, this behavior is due to the fact that the higher-order effects that we are estimating are controlled by the strong gauge coupling, and the latter decreases when the scale at which it is computed, which is proportional to the Higgs mass, increases.
Finally, an independent variation of the renormalization and factorization scales shows that, in this case, the dominant uncertainty is given by the dependence on the factorization scale.

Definition of the Higgs-bottom coupling
In the production of a SM-like Higgs boson, the contribution of bottom-quark annihilation and the effect of the bottom-quark loops in gluon fusion amount to a few percent of the total cross section. Therefore, in that case the theoretical uncertainty associated to the definition of the Higgs coupling to bottom quarks is negligible compared to other sources of uncertainty. On the other hand, this uncertainty becomes significant in scenarios where the Higgs-bottom coupling is enhanced with respect to its SM counterpart, In the MSSM the tree-level couplings of the neutral Higgs bosons to bottom quarks are modified as follows: where α is the mixing angle in the CP-even Higgs sector. In the decoupling limit, m A m Z , the mixing angle simplifies to α ≈ β − π/2, so that the coupling of h to bottom quarks is SM-like, while the couplings of H and A are both enhanced by tan β.
In this section we discuss two issues that affect the precise definition of the Higgs-bottom couplings: the first concerns the choice of renormalization scheme -and scale -for the bottom mass from which the couplings are extracted; the second concerns higher-order effects in the procedure through which the tan β-enhanced SUSY contributions are resummed in effective Higgs-bottom couplings. to a 17% decrease in the cross section for h production, and a 24% decrease in the cross section for H production. The use of m b (m φ /2) would instead decrease the cross section for h production by 34%, and the one for H production by 51%, with respect to the values obtained with M b . As a second example, we take the light-stop scenario with m A = 300 GeV and tan β = 10, where the lightest scalar h has SM-like couplings to quarks. In this case the cross section for h production varies by less than 2% when choosing among the three options discussed above for the definition of Y h b . For the heaviest scalar H, on the other hand, the changes in the cross section relative to the value derived with M b amount to −22% and −50% when Y H b is extracted from m b (m b ) and m b (m H /2), respectively. The strong sensitivity of the production of non-standard Higgs bosons on the choice of renormalization scheme (and scale) for the bottom mass and Yukawa coupling has been discussed in the past, see e.g. refs. [8,49,88]. However, unlike many other processes for which there are theoretical arguments in favor of one or the other choice, for Higgs production in gluon fusion we are not aware of any such arguments that go beyond heuristic. As was already noted in ref. [8], the options of relating Y φ b to M b or to m b (m b ) might seem preferable to the one of using m b (m φ /2), in that they lead to smaller two-loop contributions. If in the one-loop part of the amplitude for scalar production we identify the mass of the bottom quark with M b and the bottom Yukawa coupling with m b (µ b ), where µ b is a generic renormalization scale, the contribution of diagrams with bottom quarks and gluons to the two-loop part of the amplitude reads

Scheme and scale dependence of the bottom mass
where C F = 4/3 and C A = 3 are color factors, τ = 4 m 2 b /m 2 φ , and we omit an overall multiplicative factor. Truncating the functions at the first order in an expansion in powers of τ , one finds [11] with The equations above show that the two-loop bottom contribution to the gluon-fusion amplitude contains powers of ln(m 2 φ /m 2 b ), and that the choice µ b = m b does eliminate some of the logarithmically enhanced terms. Similarly, relating the coupling entering the one-loop part of the amplitude to the pole mass M b eliminates the whole piece proportional to F 1 1/2 (τ ) in eq. (5). Each of the two remaining terms, C F F C F (τ ) and C A F C A (τ ), also contains powers of ln(m 2 φ /m 2 b ), but for realistic values of m φ the two terms largely cancel out against each other, resulting in a small two-loop contribution from bottom quarks. However, such cancellation should be considered accidental: there is no argument suggesting that it persists at higher orders in QCD, or that it is motivated by some physical property of the bottom contribution to gluon fusion. To illustrate this point, we can consider the case of Higgs decay to two photons: the one-loop bottom contribution to the amplitude has the same structure as the corresponding contribution to gluon fusion, but the two-loop bottom-gluon contribution is obtained from eq. (5) by dropping the term proportional to C A , which originates from diagrams with threeand four-gluon interactions. In that case no significant cancellation occurs, and the amplitude is not In fact, it was also noted in ref. [8] that the two-loop bottom-gluon contribution to the amplitude for Higgs decay to photons is minimized when the one-loop contribution is fully expressed in terms of m b (m φ /2).
In the case of the Higgs coupling to photons, the problems related to the ambiguity in the definition of Y φ b have been solved with a resummation of the leading and next-to-leading logarithms of the ratio . Until a similar calculation is performed for the Higgs coupling to gluons, there is no obvious reason to favor one choice of renormalization scheme (and scale) for the bottom Yukawa coupling over the others. In our study we choose to relate the coupling to the pole mass M b , and we consider the difference between the results obtained using M b and those obtained using m b (m φ /2) as a measure of the uncertainty associated with the uncomputed higher-order QCD corrections. For the production of a SM-like Higgs with mass around 125.5 GeV, this procedure -also advocated by the LHC-HXSWG in ref.
[3] -results in an uncertainty of 1 -2% in the gluon-fusion cross section. On the other hand, as we show in figures 12 and 13 for scalar and pseudoscalar production in the light-stop scenario, the cross section could be reduced by more than 60% in the regions of the m A -tan β plane where the gluon-fusion process is dominated by the bottom-quark contribution. It is however worth recalling that, as shown in figures 3 and 4, in such regions the total cross section for Higgs production is dominated by bottom-quark annihilation. In the 5FS, the cross section for the latter process is known at the NNLO in QCD [43,44], and it is free of large logarithms of the ratio m 2 The theoretical uncertainty of the cross section for bottom-quark annihilation associated to reasonable variations around this scale choice is already included in the uncertainty bands shown in figures 10 and 11 in the previous section.

Resummation of tan β-enhanced corrections
It is well known that, in the MSSM, loop diagrams involving superparticles induce tan β-enhanced corrections to the couplings of the Higgs bosons to bottom quarks [90]. If all superparticles are considerably heavier than the Higgs bosons they can be integrated out of the MSSM Lagrangian, leaving behind a two-Higgs-doublet model with effective Higgs-bottom couplings where Y φ b are the tree-level Higgs-bottom couplings defined in eq. (4), and, retaining only the O(α s ) contribution from diagrams with sbottoms and gluinos, the tan β-enhanced term ∆ b reads   figure 12 for the production of the pseudoscalar A.
In the limit m A m Z , where cot α ≈ − tan β, the superparticle contributions encoded in ∆ b decouple from the coupling of the lightest scalar, while the couplings of the heaviest scalar and of the pseudoscalar are both rescaled by a factor (1 − ∆ b cot 2 β)/(1 + ∆ b ).
In refs. [39,40] it was shown that, in the calculation of processes that involve the Higgs-bottom couplings, the tan β-enhanced corrections can be resummed to all orders in the expansion in powers of ∆ b by inserting the effective couplings of eq. (10) in the lowest-order amplitude for the considered process. In the case of gluon fusion, this amounts to using Y b φ in the bottom contribution to the one-loop part of the amplitude. However, when this resummation procedure is combined with the actual calculation of the superparticle contributions to the one-and two-loop amplitude for gluon fusion, care must be taken to avoid double counting. To this effect, we must subtract from the full result for the two-loop amplitude the contribution obtained by replacing Y b φ in the resummed one-loop amplitude with the O(∆ b ) term of the expansion of Y b φ in powers of ∆ b . Depending on the choice of renormalization scheme for the parameters in the sbottom sector, additional tan β-enhanced terms could be induced in the two-loop amplitude by the counterterm of the Higgs-sbottom coupling that enters the sbottom contribution to the one-loop amplitude. To avoid the occurrence of large two-loop corrections, which would put the validity of the perturbative expansion into question, we employ for the sbottom sector the OS renormalization scheme described in ref. [32].
An ambiguity in the procedure for the resummation of the ∆ b terms concerns the treatment of the Higgs-bottom couplings entering the two-loop part of the gluon-fusion amplitude. The difference between the results obtained using either Y φ b or Y b φ in the two-loop part is formally of higher order, i.e., it amounts to three-loop terms that are suppressed by a factor A b /(µ tan β) with respect to the dominant three-loop terms of O(∆ 2 b ) accounted for by the resummation. Nevertheless, in our study we choose to identify the Higgs-bottom couplings in both the one-and two-loop parts of the amplitude with Y b φ . We found that this choice allows us to reproduce -after an expansion in powers of ∆ bthe three-loop result that can be inferred from ref. [40], where the sub-dominant terms proportional to A b were also resummed in the effective couplings.
For large values of tan β, the factor ∆ b can even become of order one, unless the superpotential parameter µ is suppressed with respect to the soft SUSY-breaking masses. The effect of the SUSY correction on the effective Higgs-bottom couplings depends crucially on the sign of ∆ b . For positive ∆ b the correction suppresses the couplings, reducing the overall relevance of the bottom contribution to gluon fusion. On the other hand, for negative ∆ b the correction enhances the couplings, which diverge as ∆ b approaches −1. As a consequence, when ∆ b is large and negative the result for the gluon-fusion cross section is extremely sensitive to the precise value of ∆ b , and a refined calculation of the latter becomes mandatory to reduce the uncertainty associated to the bottom contribution.
The first obvious step to improve the calculation of ∆ b consists in including other one-loop contributions that are not shown in eq. (11). In particular, the diagrams with stops and charginos induce a contribution, controlled by the top Yukawa coupling, that can be comparable in size with the O(α s ) contribution in eq. (11). In our numerical analysis we use by default the full one-loop result for ∆ b as computed by FeynHiggs, which allows us to resum in our prediction for the Higgs-production cross section also the tan β-enhanced corrections of electroweak origin.
Another improvement in the calculation would come from the inclusion of the dominant two-loop contributions to ∆ b , which have been computed in ref. [41] but are not yet implemented in FeynHiggs. Indeed, it was shown in ref. [41] that the one-loop result for ∆ b is particularly sensitive to changes in the renormalization scales at which the strong-gauge and top-Yukawa couplings are expressed, and that the inclusion of the two-loop contributions stabilizes this scale dependence. In particular, both the one-loop sbottom-gluino and stop-chargino contributions to ∆ b vary by roughly ±10% when the renormalization scales are lowered or raised by a factor of two around their central values, which are chosen as the average of the masses of the relevant superparticles. We can therefore estimate the uncertainty of the gluon-fusion cross section associated to the one-loop computation of ∆ b by varying by ±10% the result provided by FeynHiggs.
In general, the impact of the uncertainty of ∆ b on the total uncertainty of the gluon-fusion cross section depends on the considered point in the MSSM parameter space. As was the case also for the scheme and scale dependence of Y φ b discussed in the previous section, the ∆ b uncertainty can be significant only if the bottom contribution to the cross section is substantially enhanced with respect to the SM case. For illustration, we consider again the light-stop scenario with m A = 130 GeV and tan β = 40, where both Higgs scalars have enhanced couplings to bottom quarks. The superpotential parameter µ has positive sign, and the ∆ b corrections suppress the effective couplings Y b φ . We find that the cross sections for h and H production in gluon fusion increase by 4% and 7%, respectively, if the value of ∆ b is reduced by 10%, while they decrease by 4% and 6%, respectively, if ∆ b is increased by 10%. The effect is larger if µ is taken negative, so that the ∆ b corrections enhance the effective couplings. In that case the dependence on ∆ b is reversed: if we consider the same point in the lightstop scenario but flip the sign of µ, the cross sections for h and H production in gluon fusion decrease by 17% and 16%, respectively, when |∆ b | is reduced by 10%, while they increase by 23% and 21%, respectively, when |∆ b | is increased by 10%.
Finally, we stress that a similar uncertainty affects the cross section for Higgs production via bottom-quark annihilation, where the tree-level amplitude is computed in terms of the effective couplings Y b φ . Also in this case, we can estimate the uncertainty by varying by ±10% the value of ∆ b provided by FeynHiggs.

Uncertainties from the PDFs and α s
The prediction for the total cross section at hadron level is affected by our imperfect knowledge of the proton PDFs. This uncertainty has different sources: the PDFs cannot be computed from first principles but they rather have to be fitted from data, and the experimental error of the latter affects the outcome of the fit and propagates to the prediction of any observable. Also, the choices related to the fitting methodology and to the mathematical representation of the PDFs induce an ambiguity in the results, as can be appreciated by comparing the PDF parameterizations provided by three collaborations that perform a global fit of low-and high-energy data: MSTW2008 [83], CT10 [91] and NNPDF2.3 [92]. These uncertainties will be discussed in section 3.3.1, together with the parametric dependence of the cross section on the value of the strong coupling constant. Another source of uncertainty is related to the available perturbative-QCD information on the scattering processes from which the PDFs are extracted. Among these perturbative effects, an issue that is particularly relevant in the case of Higgs production via bottom-quark annihilation is the consistent inclusion of the bottommass effects in the evolution of the PDFs according to the DGLAP equations. The transition between four and five active flavors in the proton occurs at a matching scale that is set equal to the bottom mass. The bottom density in the proton depends parametrically on this matching scale, which in turn affects the predictions for the cross section. The phenomenological implications of this issue will be discussed in detail in section 3.3.2. A systematic discussion of further sources of theoretical uncertainty -such as, e.g., the dependence of the PDFs on the choice of renormalization and factorization scale in the matrix elements that are used to perform the fit -is not yet available in the literature, and goes beyond the scope of this paper.

Combination of PDF and α s uncertainties
The uncertainty associated to the experimental errors of the data from which the PDFs are extracted is represented by the PDF collaborations with the introduction of N R different PDF sets (replicas), all equivalent from the statistical point of view in the description of the data. Any observable has to be computed N R times with the different sets, and the spread of the results can be interpreted as the error induced by the PDF due to the data and to the fitting methodology. The replicas are determined by the PDF collaborations following the Hessian (for MSTW2008 and CT10) or the Monte Carlo (for NNPDF2.3) approaches, and the PDF error has to be computed accordingly. In QCD the cross sections are also affected by a parametric uncertainty associated to the input value of the strong coupling constant. This dependence is particularly relevant in the gluon-fusion cross section, which is proportional to α 2 s at the LO and is subject to very large QCD corrections, of O(α 3 s ), at the NLO. Each PDF collaboration recommends a different central value for α s (m Z ), generating a spread of the central predictions for the Higgs-production cross section. The combination of the PDF and α s uncertainties (henceforth, PDF+α s ) and their correlation was first discussed in ref. [93]. A conservative approach to combine the different predictions obtained using the MSTW2008, CT10 and NNPDF2.3 PDF sets is known as PDF4LHC recipe, and it amounts to taking the envelope of the PDF+α s uncertainty bands of the three collaborations, where for each group the preferred α s (m Z ) central value is adopted [94]. Following this reference we take ∆α s = ±0.0012 for the experimental error on the strong coupling constant.
Due to the very steep behavior of the PDFs for increasing values of the final-state invariant mass, the gluon-fusion process receives its dominant contribution from the threshold production region, with a very important role played by the virtual corrections and by the universal, factorizable, soft-gluon corrections. Consequently, the cross section is dominated by the LO-kinematics configurations also at higher perturbative orders. At the LO, the gluon-fusion cross section depends on the rapidity of the Higgs boson only through the PDFs, therefore the relative size of the PDF+α s uncertainty does not depend on the details of the partonic process, but only on the value of the Higgs-boson mass. As a consequence, the relative PDF+α s uncertainty, for a given value of the Higgs mass, can be read directly from the tables of the SM predictions reported in the appendix B of the latest LHC-HXSWG report [5]. Differences with respect to the SM predictions may originate from hard, process-dependent radiative corrections, but their impact on the relative PDF+α s uncertainty is at the sub-percent level.
To assess the PDF+α s uncertainty of the cross section for Higgs production in bottom-quark annihilation we adopt again the PDF4LHC recipe. The bottom density in the proton does not have an intrinsic component, but it is generated dynamically, via gluon splittings, by the DGLAP evolution of the PDFs. Therefore, the uncertainties of the bottom and gluon PDFs are strongly correlated.
Similarly to the case of gluon fusion, for a given value of the Higgs mass the relative PDF+α s uncertainty of the cross section for bottom-quark annihilation differs very little between the SM and the MSSM, because the radiative corrections involving SUSY particles affect the kinematics of the process only at higher orders. 6 We find that the uncertainty has an almost constant behavior when the mass m φ of the produced Higgs boson is lighter than 300 GeV, and that it increases for larger mass values: for example, at the NNLO, the PDF+α s uncertainty of the cross section for bottom-quark annihilation amounts to ± 6/6/8/21% for m φ = 124/300/500/1000 GeV.

Bottom-mass dependence of the PDFs
The calculation of hadronic cross sections involves the convolution of the partonic cross sections with the PDFs, which have an intrinsic dependence on the bottom mass. For example, the central set of MSTW2008 [83], which we use as default for our analysis, assumes a pole mass M b = 4.75 GeV. Converted to the MS mass via a three-loop QCD calculation, this corresponds to m b (m b ) = 4.00 GeV, which differs both from the value recommended by the LHC-HXSWG, m b (m b ) = 4.16 GeV [57,59], and from the current PDG value, m b (m b ) = 4.18 GeV [95].
In addition to their dependence through the PDFs, the cross sections for Higgs production also depend on the bottom mass at the partonic level, i.e., through the bottom Yukawa coupling, the bottom-quark propagators and the phase space. In the regions of the MSSM parameter space where the bottom-quark contributions to Higgs production are enhanced, it becomes vital to evaluate the partonic cross sections with the correct input value for the bottom mass, which, as mentioned above, may not necessarily correspond to the value used in the PDFs. In this section we will examine the uncertainty that arises when we choose the bottom mass entering the partonic cross sections independently from the PDF set.
The MSTW2008 PDFs come in seven sets obtained with M b ranging from 4 GeV to 5.5 GeV in steps of 0.25 GeV. In ref. [96] the MSTW collaboration studied the sensitivity of the PDFs on the value of the bottom mass, showing that the PDFs for the gluon and for the four lightest quarks are almost insensitive to M b , whereas the bottom PDF exhibits quite a strong dependence. As shown in figure 6 of ref. [96], a variation by ±0.5 GeV around the central value M b = 4.75 GeV leads to changes in the bottom PDF that exceed the 90% C.L. uncertainty, even for the relatively large value of the factorization scale relevant to Higgs production, µ F ≈ 100 GeV.
The cross section for Higgs production via gluon fusion is mostly sensitive to the gluon PDF, and receives only a small contribution, starting at the NLO, from diagrams with initial-state bottom quarks. As a result, when we evaluate the gluon-fusion cross section with the seven PDF sets -while fixing the bottom mass in the partonic cross section -we find that the result changes only at the per mil level, independently of the phenomenological scenario under consideration. We conclude that, for this process, the formal inconsistency of choosing different values for the bottom mass in the partonic cross section and in the PDFs induces only a negligible uncertainty.
In contrast, the hadronic cross section for Higgs production in bottom-quark annihilation, when computed in the 5FS, depends directly on the bottom PDF. As a result, we expect this process to show a significant dependence on the value of the bottom mass used in the PDFs, and the issue of consistency with the bottom mass used in the definition of the bottom Yukawa coupling becomes unavoidable.
In figure 14 we investigate the bottom-mass dependence of the hadronic cross section for pseudoscalar production in bottom-quark annihilation (we find similar behaviors for the production of the scalars, both light and heavy). The plot on the left shows the hadronic cross section as a function of the pseudoscalar mass m A , in the light-stop scenario with tan β = 20. As in section 2, the renormalization and factorization scales are set to µ R = m A and µ F = m A /4. The central (black) solid line in the left plot is computed with our default settings, namely we use the PDF set with M b = 4.75 GeV and we relate the Yukawa coupling Y b A to m b (m A ), which we obtain from the input m b (m b ) = 4.16 GeV via renormalization-group evolution. The plot on the right of figure 14 represents the variation of the cross section relative to this default setting, when we change the bottom mass in the PDFs and/or in the Yukawa coupling.
In both plots, the red band between dot-dashed lines indicates the spread in the cross section obtained with the extreme PDF sets -corresponding to M b = 4 GeV and M b = 5.5 GeV, respectively -with Y b A fixed to the default value. As expected, the impact of the bottom mass used in the PDFs is significant: it amounts to about (+20/-15)% at large m A , with larger values of M b corresponding to smaller cross sections. This anti-correlation is a consequence of the fact that, for larger bottom masses, the reduced available phase space for the splitting of gluons into bottom pairs leads to a suppression of the bottom PDF. On the other hand, the bottom Yukawa coupling is directly correlated with the magnitude of the cross section. Simultaneously adjusting the bottom mass entering the bottom Yukawa coupling and the one entering the bottom PDF should therefore lead to some degree of compensation between these two effects.
Converting , shown as a blue band between solid lines in both the left and the right plots, is thus also asymmetric, and it is of the order of (+15/-30)% at large m A .
Our procedure to estimate the uncertainty of the cross section for bottom-quark annihilation arising from the bottom-mass dependence of the PDFs is similar to the one in ref. [88]. We fix the bottom Yukawa coupling to the value implied by m b (m b ) = 4.16 GeV, as recommended by the LHC-HXSWG, and we use as uncertainty the spread in the cross section caused by the variation of M b in the PDFs around the central value of 4.75 GeV. However, the full variation of ±0.75 GeV allowed by the MSTW2008 PDFs, which would correspond to the red band in figure 14, seems overly conservative for our purposes. A variation of ±0.25 GeV is in fact sufficient to encompass the value M b = 4.92 GeV, which corresponds at the three-loop level to the recommended MS mass m b (m b ) = 4. 16 GeV. This variation finally leads to an estimate of the uncertainty of about ±6%. A similar estimate is obtained from NNPDF2.1 [97], which also provides PDF sets with different values of M b .

Higher-order SUSY contributions to gluon fusion
In this section we discuss two sources of uncertainty affecting the SUSY contributions to the cross section for gluon fusion. The first is the validity of the expansion in the heavy superparticle masses of the two-loop SUSY contributions; the second is the impact of the three-loop SUSY contributions that are not included in SusHi.

Validity of the expansion in the SUSY masses
The results implemented in SusHi for the two-loop stop contributions to lightest-scalar production rely on the VHML, while the results for the remaining two-loop SUSY contributions rely on expansions in inverse powers of the superparticle masses. The latter include terms up to where m φ denotes a Higgs mass and M denotes a generic superparticle mass. Therefore, the validity of the results for the two-loop SUSY contributions is limited to the region where the mass of the produced Higgs boson is smaller than the lowest-lying SUSY-particle threshold of the Feynman diagrams involved. In all of the six benchmark scenarios considered in our study, the lightest-scalar mass lies comfortably below this limit. Since we consider m A ≤ 1 TeV, the same applies also to the masses of the heaviest scalar and of the pseudoscalar in the five scenarios in which the squark masses are themselves of the order of 1 TeV. In the light-stop scenario, on the other hand, the lowest-lying SUSY threshold is at 2 mt 1 ≈ 650 GeV, hence our need to limit our analysis to m A ≤ 500 GeV.
To assess the quality of our approximation in the vicinity of the threshold, we multiply the twoloop stop and sbottom contributions to the gluon-fusion amplitude by test factors tq ≡ A 1 q 1 /A 1 , exp q 1 , withq = {t,b}. Specifically, A 1 q 1 is the lightest-squark contribution to the one-loop part of the scalarproduction amplitude including the full mass dependence, while A 1 , exp q 1 includes only the leading O(m −2 q 1 ) terms in the expansion in the lightest-squark mass. Assuming that the expanded two-loop contributions deviate from the full ones by an amount comparable to that seen in the one-loop contributions, the variation in the gluon-fusion cross section resulting from the introduction of the test factors can be considered as an estimate of the uncertainty associated to the expansion in the SUSY masses.
The contour plots in figure 15 show the effect of introducing these test factors on the cross section for the production of the heaviest scalar (left plot) and of the pseudoscalar (right plot) in the light-stop scenario. In the case of H production, the variation of the cross section at large m A amounts to a few percent when tan β is sufficiently large, but it can exceed 20% when 8 < ∼ tan β < ∼ 16. As can be seen in the right plot of figure 5, in this region the one-loop quark and squark contributions to the gluonfusion amplitude largely cancel each other, with the result that the total cross section becomes small and particularly sensitive to variations in the two-loop contributions. This sensitivity to higher-order effects manifests also as the large scale uncertainty, up to 50%, visible in the right plot of figure 8. In the case of A production, on the other hand, no such cancellations occur, because the squarks do not contribute to the one-loop amplitude for gluon fusion. 7 The variation of the cross section at large m A is therefore limited to a few percent even for moderate tan β. We performed the same analysis on the other five benchmark scenarios, where the squark masses are of the order of 1 TeV. As expected, we found that the effect of rescaling the two-loop SUSY contributions by test factors tq is much smaller than in the light-stop scenario, and it is certainly negligible when compared to the scale uncertainty of the cross section. In particular, in the tau-phobic scenario -where the squark contributions to the gluon-fusion amplitude are enhanced by the large value of the parameter µ -the effect on H production reaches the few-percent level only when m A approaches 1 TeV, for moderate tan β. In the remaining four scenarios the effect is even smaller.

The SUSY contributions at the NNLO
The QCD corrections to the gluon-fusion cross section are large, typically exceeding 100% at the energy of the LHC. In the SM, an excellent approximation to these corrections is obtained in the VHML (or heavy-top limit) [12,13,14], where a perturbative K-factor is calculated in the effective theory that results from neglecting the bottom Yukawa coupling and integrating out the top quark, leaving behind a point-like Higgs-gluon interaction term L ggH = −(1/4v) C(α s ) HG µν G µν , with v ≈ 246 GeV. The Wilson coefficient accounts for heavy particles that mediate the Higgs-gluon coupling in the underlying theory. In the SM, this is just the top quark; it is easy to see, though, that the inclusion of stop squarks (and gluinos) only affects C(α s ), while the form of L ggH remains unchanged. A comparison with the full result at the NLO suggests that, within the SM, the VHML provides a decent approximation of the NNLO top contributions also for rather large Higgs masses [81,82]. Therefore, SusHi includes the NNLO effects in the cross sections for the production of all three neutral Higgs bosons of the MSSM. Within the effective theory, the K-factor at the NNLO takes the form where Σ (n) is the n th -order term in the perturbative expansion of the hadronic cross section based on L ggH C(αs)≡1 . Note that, in the NNLO part of the K-factor in eq. (13), the only genuine three-loop term that depends on the underlying theory is C (2) . This observation was exploited in ref. [36] to derive an estimate of the NNLO top/stop contribution to the gluon-fusion cross section in the MSSM. In particular, it was shown that the final result depends only very weakly on the numerical value of C (2) . Consequently, once the two-loop stop contributions are included in C (1) , the unknown three-loop stop contributions to C (2) induce an uncertainty in the cross section much smaller than the residual uncertainty derived from scale variation. It was suggested to use the top contribution C (2) t for the whole C (2) , and to estimate the related uncertainty by varying that coefficient within the interval [0, 2 C In ref. [36] the hadronic cross section was obtained, in analogy to the SM NNLO result, by reweighting its exact LO expression with the K-factor of eq. (13): where A 1 tt ≡ A 1 t + A 1 t is the one-loop amplitude including both the top and stop contributions with the exact Higgs-mass dependence (in particular, A 1 tt → C (0) in the VHML, i.e. for m φ → 0). However, as was discussed also in the previous section, there exist so-called gluophobic regions of the MSSM parameter space in which the top and stop contributions to the amplitude can cancel each other to a large extent. Since the precise values of the MSSM parameters where this cancellation is maximal differ between the full calculation and the VHML, the ratio |A 1 tt |/C (0) entering the cross sectionsee eqs. (13) and (14) -can become spuriously large when C (0) ≈ 0. In order to evade this effect, we replace C (0) in eq. (12) with A 1 tt . This leads to the following expression for the cross section: This formula applies to both MSSM scalars. The effective Lagrangian for the gluonic interaction of the pseudoscalar involves an additional operator which contributes at the NNLO [35], but it can be treated in a completely analogous way.
In SusHi, the NNLO top and stop contributions to the gluon-fusion cross section in the VHML are isolated by subtracting from the σ NNLO in eq. (15) the same quantity truncated at the NLO (and computed with NLO PDFs). The result is then added to the full NLO cross section, which accounts also for the bottom and sbottom contributions and for the known Higgs-mass dependence of the two-loop amplitude. The 6% suppression of the cross section for the production of a SM-like scalar induced by the NNLO stop contributions in the light-stop scenario -see section 2.2 -can be ascribed to the effect of the term 2 C (1) Σ (1) Re A 1 tt in the second line of eq. (15). Indeed, the large value of the (normalized) NLO term of the cross section in the effective theory, Σ (1) /Σ (0) ≈ 26, compensates for the suppression by α s /π, with the result that the effect of the two-loop stop contribution to C (1) at the NNLO is roughly as large as the corresponding effect at the NLO.
To assess the uncertainty arising from the fact that we neglect the three-loop SUSY contributions to C (2) , we make use of a recent calculation of those contributions in the VHML [37,38]. The calculation is based on an expansion of the relevant Feynman diagrams in terms of certain hierarchies among the different masses, similar to the strategy that was pursued in ref. [70] for the calculation of the 3-loop corrections to the Higgs mass in the MSSM. The results of ref. [38] are available in the form of a Mathematica file, which provides the basis for the expansion of C (2) in various hierarchies of the masses mt 1 , mt 2 , mg, m t , and mq, combined with expansions in differences of these masses. Following an algorithm suggested in ref. [38], these expansions should allow one to derive a numerical approximation for C (2) in any viable MSSM scenario.
Applying this approach to the scenarios defined in table 1, we find that the deviation of the whole C (2) from the top contribution C (2) t is rather small, and the second-order coefficient certainly stays within the range [ 0, 2 C (2) t ]. Varying C (2) within this interval, we estimate that the effect of the threeloop SUSY contributions to the gluon-fusion cross section does not exceed 1% in all of the scenarios considered in this paper. It is therefore a viable strategy to follow ref. [36] and set C (2) = C (2) t , attributing an uncertainty of ±1% to the final result for the cross section.

Conclusions
A precise prediction of the cross sections for Higgs-boson production, as well as a detailed understanding of the associated uncertainties, are of vital importance to interpret the recent discovery of a Higgs boson in the context of the MSSM. In this paper we used the public code SusHi [53], which computes the cross sections for gluon fusion and bottom-quark annihilation, to study the production of scalar and pseudoscalar Higgs bosons in a set of MSSM scenarios compatible with the LHC results. We showed how the cross sections can substantially differ from the SM prediction, and how their qualitative behavior over the MSSM parameter space depends mainly on the relative importance of the contributions involving top and bottom quarks. We also emphasized that, in a scenario with relatively light squarks which is not yet constrained by the LHC, the contributions to the gluon-fusion process that involve superparticles can significantly suppress the cross section for scalar production.
Next, we studied the different sources of uncertainty that affect our predictions for the Higgsproduction cross sections. Some of these uncertainties, namely the ones associated to the choice of renormalization and factorization scales, to the PDF parameterization and to the input value for the strong coupling constant, are relevant also for the production of the SM Higgs, although their size may differ in the case of the production of non-standard Higgs bosons. In contrast, the uncertainties associated to the definition of the bottom mass and Yukawa coupling are practically negligible in the SM -where the bottom-quark contributions amount only to a few percent of the total cross sectionbut they can become dominant in regions of the MSSM parameter space where the couplings of the Higgs bosons to bottom quarks are enhanced. In the particular case of heavy-scalar and pseudoscalar production at large tan β, we found that legitimate variations in the renormalization scheme and scale of the bottom Yukawa coupling can suppress the gluon-fusion cross section by more than 60%, due to the presence of large QCD corrections enhanced by logarithms of the ratio m 2 φ /m 2 b . Luckily, in this case the total cross section is dominated by the contribution of bottom-quark annihilation, which is subject to a considerably smaller scale uncertainty. Finally, we studied the uncertainties associated to our implementation of the SUSY contributions to gluon fusion at the NLO and, partially, at the NNLO. With the exception of a gluophobic region in the light-stop scenario, these uncertainties are generally small, reflecting the sub-dominant nature of the SUSY contributions themselves for values of the squark masses compatible with the LHC bounds.
Future improvements in the accuracy of the predictions for the Higgs-production cross sections in the MSSM could come from different directions. First of all, any progress in the SM calculation will eventually trickle down to the MSSM calculation, at least where the production of a SM-like scalar is concerned. In addition, a resummation of the QCD corrections enhanced by ln(m 2 φ /m 2 b ), analogous to the one performed in ref. [89] for the Higgs decay to photons, will be necessary to reduce the large uncertainty in the production of non-standard Higgs bosons via gluon fusion (incidentally, such calculation would benefit all models with enhanced Higgs couplings to bottom quarks, whether they are supersymmetric or not). Implementing the existing results for the two-loop contributions to ∆ b [41], in both the Higgs mass and cross-section calculations, will also reduce the uncertainty in scenarios where the bottom contributions are relevant. Finally, it could be worthwhile to improve the calculation of the gluon-fusion cross section by taking into account the full Higgs-mass dependence of the two-loop squark-gluon 8 contributions [10,11,24] -to cover scenarios in which the non-standard Higgs bosons are heavier than the third-generation squarks -and by including the genuine three-loop effects [37,38].
For the time being, however, we believe that SusHi provides the most sophisticated calculation of the Higgs-production cross sections in the MSSM available to the physics community. Differently from the case of the SM, where all the relevant inputs are now known and a definite prediction for the total cross section can be made solely as a function of the collision energy, in the MSSM the predictions for the cross sections -and the relevance of the different sources of uncertainty -depend crucially on a number of yet-undetermined parameters. In the appendix we collect predictions for the Higgs-production cross sections via gluon fusion and bottom-quark annihilation, and their respective uncertainties, in a few representative points of the six benchmark scenarios described in section 2.1. However, the tables in the appendix should be regarded as having an illustrative purpose only: we encourage the readers to take both SusHi and our recipes for the uncertainties directly in their own hands, and use them to analyze their favorite corners of the MSSM parameter space. Our results should prove useful for ruling out scenarios that are incompatible with the current experimental bounds. We also hope that, when the time comes, they will help interpret within the MSSM the discovery of new particles at the LHC.

Acknowledgments
This work was initiated in the context of the activities of the LHC-HXSWG, and some of our preliminary findings were presented in section 14 of ref.
[5]. We thank the authors of ref. [56] and the SUSY convenors of the ATLAS and CMS collaborations for helpful communication about the experimental constraints on the MSSM benchmark scenarios, and the authors of ref. [38] for helpful communication about their results. We also thank M. Spira for useful discussions about the uncertainties of the cross-section calculations. This

Appendix: Cross sections and uncertainties
In this appendix we include eighteen tables, listing the cross sections and uncertainties for the production at the LHC of the three neutral Higgs bosons in the six MSSM scenarios defined in table 1. We use version 1.3.0 of SusHi, and provide separate results for gluon fusion and bottom-quark annihilation. Input files for the six scenarios can be found on the code's website [53]. We set • The uncertainty δ∆ b , stemming from the resummation of tan β-enhanced corrections to Y φ b and described in section 3.2.2, is computed by adding an uncertainty of ±10% to the value of ∆ b obtained from FeynHiggs.
The PDF uncertainties are not included in the tables, but they were extensively discussed in section 3.3. In section 3.3.1 we pointed out that the relative size of the PDF+α s uncertainty depends mainly on the value of the Higgs mass, thus it can be taken over directly from the existing estimates for the production of the SM Higgs. Apart from the PDF+α s uncertainty, in the case of bottomquark annihilation an additional uncertainty of ±6% has to be added due to the dependence of the bottom-quark PDF on the pole bottom mass (see section 3.3.2). The uncertainties associated to our incomplete knowledge of the SUSY contributions to the gluon-fusion cross section can become sizeable only in the light-stop scenario, especially in the case of H production at large m A and moderate tan β. We do not include them in the tables, pointing the reader to the discussion in section 3.4.
We show the results for the m max