1 Introduction

CERN-PH-TH/2014-048,

DESY 14-040, IFUM 1024-FT,

LPN14-061, MCnet-14-07,

WUB/14-01

Towards precise predictions

[7pt] for Higgs-boson production in the MSSM
E. Bagnaschi, R.V. Harlander, S. Liebler, H. Mantler, P. Slavich

[7pt] and A. Vicini

LPTHE, UPMC Univ. Paris 06, Sorbonne Universités, 4 Place Jussieu, F-75252 Paris, France

LPTHE, CNRS, 4 Place Jussieu, F-75252 Paris, France

Fachbereich C, Bergische Universität Wuppertal, Gaußstraße 20, D-42119 Wuppertal, Germany

Universität Hamburg, Luruper Chaussee 149, D-22761 Hamburg, Germany

CERN, Theory Division, CH-1211 Geneva 23, Switzerland

Dipartimento di Fisica, Università di Milano and INFN, Sezione di Milano,

Via Celoria 16, I–20133 Milano, Italy

1 2 3

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.

## 1 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-to-leading 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 next-to-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 doublets, and , whose relative contribution to electroweak symmetry breaking is determined by the ratio of vacuum expectation values of their neutral components, . The spectrum of physical Higgs bosons is richer than in the SM, consisting of two neutral scalars, and , one neutral pseudoscalar, , and two charged scalars, . 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 . 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 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 top-stop-gluino contributions to the production of the lightest scalar , it might not be applicable to the production of the heaviest scalar and of the pseudoscalar , 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 -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 and to bottom quarks are enhanced by with respect to the corresponding coupling of the SM Higgs, while their couplings to top quarks are suppressed by . When that is the case, the bottom-quark contributions to the gluon-fusion process – which for a SM-like Higgs with mass around  GeV amount to roughly 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 , 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 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 -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 -enhanced squark contributions in an effective Higgs-bottom 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 Higgs-boson 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 -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.

## 2 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.

In view of these considerations, we will focus on the set of MSSM scenarios compatible with the LHC findings that has recently been proposed in ref. [56]. We will study the effect of the different contributions to the total cross section for the production of the MSSM Higgs bosons, relying on the approximate NNLO-QCD calculations implemented in SusHi.

### 2.1 The benchmark scenarios

The SM parameters entering our calculations include the -boson mass  GeV, the -boson mass  GeV, the Fermi constant and the strong coupling constant .4 For the masses of the top and bottom quarks we take the pole mass  GeV [58] and the SM running mass (in the scheme)  GeV [59].

At the tree level, the MSSM neutral scalar masses and and the scalar mixing angle can be computed in terms of , and the pseudoscalar mass only. However, the radiative corrections to the tree-level predictions can be substantial, and they bring along a dependence on all of the other MSSM parameters. To compute the masses and the couplings of Higgs bosons and superparticles in a given point of the MSSM parameter space we use the public code FeynHiggs [60], which includes the full one-loop [61] and dominant two-loop [62, 63, 64, 65, 66] corrections to the neutral Higgs masses. Since the theoretical uncertainty of the Higgs-mass calculation in FeynHiggs has been estimated to be of the order of 3 GeV [67],5 we consider as phenomenologically acceptable the points in the MSSM parameter space where FeynHiggs predicts the existence of a scalar with mass between  GeV and  GeV and with approximately SM-like couplings to gauge bosons.

In addition to and , the MSSM parameters most relevant to the prediction of the masses and production cross sections of the Higgs bosons are: the soft SUSY-breaking masses for the stop and sbottom squarks, which for simplicity we set all equal to a common mass parameter ; the soft SUSY-breaking gluino mass ; the soft SUSY-breaking Higgs-squark-squark couplings and ; the superpotential Higgs-mass parameter . In our convention for the sign of the latter, the left-right mixing terms in the stop and sbottom mass matrices are and , respectively. It should be noted that in our analysis the soft SUSY-breaking squark masses and trilinear couplings are expressed in an “on-shell” (OS) renormalization scheme, as described in refs. [62, 63] for the stop sector and in refs. [64, 65, 32] for the sbottom sector. Since the two-loop calculation of the Higgs masses implemented in FeynHiggs and the NLO-QCD calculation of the production cross section implemented in SusHi employ the same OS scheme, the input values of the soft SUSY-breaking parameters can be passed seamlessly from the Higgs-mass calculation to the cross-section calculation. Concerning the parameters , and , their definition is relevant to the Higgs-mass calculation only. In particular, and are expressed in the scheme, at a renormalization scale that FeynHiggs takes by default equal to , while is identified with the pole mass of the pseudoscalar. Finally, the choice of renormalization scheme for 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 , 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 , , and for the soft SUSY-breaking wino mass are listed in table 1. We vary the parameters and within the ranges

 2≤tanβ≤50 ,             90 GeV≤mA≤1 TeV . (1)

In all scenarios the Higgs-sbottom-sbottom coupling is set equal to , the left-right mixing of the first-two-generation squarks is neglected and the bino mass is obtained from the GUT relation , with the exception of the fourth scenario where we set  GeV.6 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  GeV and  GeV; the sbottom masses depend on , but the lightest sbottom is always heavier than  GeV, while the heaviest one is always lighter than  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 stop- and sbottom-pair production, based on the full 8-TeV data sample, considering the decay chains

[76, 77] ,             [76, 77] ,             [78, 79] ,

[73, 80] ,                            [71, 74] .

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, GeV, together with GeV, the masses for the lightest chargino and neutralino have a mild dependence on , but they stay within the ranges GeV and GeV for . In this case the lightest stop decays almost entirely through the loop-induced, flavor-violating channel . This channel has been investigated by ATLAS [78] and CMS [79], but the resulting bounds only reach to values of around  GeV. For the lightest sbottom, the two-body decays and (with up to 3 or 4) are kinematically open. The direct decay of to the lightest neutralino would be constrained by the searches in refs. [73, 80], but that channel is never dominant in the considered range of parameters and the experimental bounds only reach to values of below GeV. Finally, the heaviest stop and sbottom can decay through a multitude of channels, and their direct decays to or are significantly suppressed.

### 2.2 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,7 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, , 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  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  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   [13, 84], where denotes the considered Higgs boson. For bottom-quark annihilation, the central values of the scales are chosen as and   [43, 44, 85]. In the calculation of the gluon-fusion cross section we relate the bottom Yukawa coupling to the pole mass , computed at the three-loop level [86] from the input value for the running mass, . In the case of bottom-quark annihilation, on the other hand, we relate the bottom Yukawa coupling to , in turn obtained from via four-loop renormalization-group evolution [87]. In both cases, the -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 () and of the pseudoscalar (), respectively, as contour plots in the 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 (left) and (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  GeV at large . The heaviest-scalar mass grows with , and we show only the contour corresponding to  GeV to avoid clutter (for large , the contours are roughly at and independent of ). The -axis of the plot for production ends at  GeV because, for larger values, the cross section becomes essentially independent of . The -axis of the plots for and ends at  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 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, , the lightest scalar has SM-like couplings to quarks, while its mass is essentially independent of and, for , depends only weakly on . The cross section for 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  GeV, on the other hand, the couplings of to top (bottom) quarks are non-standard, being suppressed (enhanced) by . In this narrow region the total cross section for production is dominated by the contributions of the diagrams that involve the Higgs-bottom coupling, and it grows significantly with .

The behavior of the cross section for production in the plane (right plot in figure 1) is different from – and somewhat complementary to – the one for production. In the strip where  GeV, the heaviest scalar has a mass around GeV and significant couplings to both top and bottom quarks, and the cross section for its production grows with . For larger , on the other hand, grows together with , and the couplings of to top (bottom) quarks are suppressed (enhanced) by . The total cross section for production is therefore dominated, already for moderate , by the contributions of the diagrams that involve the Higgs-bottom coupling. The latter grow significantly with , but decrease with , being suppressed by powers of the ratio . Finally, the pseudoscalar couplings to top (bottom) quarks are suppressed (enhanced) by for all values of . Therefore, the behavior of the cross section for production in the plane, see figure 2, resembles the behavior of production when  GeV, and the one of production for larger : in both cases, the cross section grows with , but decreases with .

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 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 is large enough that the couplings of the lightest scalar are SM-like, gluon fusion is by far the dominant process for production, and the contribution of bottom-quark annihilation amounts only to a few percent. Only in the strip with  GeV and , where the coupling of to bottom quarks is sufficiently enhanced by , does bottom-quark annihilation become the dominant process. Conversely, bottom-quark annihilation gives the largest contribution to the cross section for production (right plot in figure 3) when  GeV and , while in the case of production (figure 4) the cross section is dominated by bottom-quark annihilation already for  GeV, as long as .

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 Higgs-quark 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 production by as much as in the decoupling region with large and . Remarkably, in this region the partial NNLO-QCD contributions from stop loops that we include following ref. [36] account by themselves for a suppression of the cross section. The theoretical uncertainty associated to these contributions will be discussed in section 3.4. For what concerns production (right plot of figure 5), the squark contributions reduce the cross section by up to for low values of , and the suppression becomes even stronger with increasing pseudoscalar mass. In particular, near the lower-right corner of the plot, where  GeV and ranges between and , the interference between the quark and squark contributions induce a suppression of the cross section by . In this region the top contribution is suppressed by , while the bottom contribution is suppressed by and only moderately enhanced by , 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 production, the effect of the squark contributions on the cross section for gluon fusion in the light-stop scenario is always less than . 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 Higgs-gauge boson couplings.8 In figure 7 we investigate the impact of the light-quark EW contributions on the production of the scalars and , plotting the ratio of the gluon-fusion cross sections computed with and without those contributions, in the 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 production (left plot) the EW corrections become fairly constant, around , in the region of sufficiently large where the lightest scalar has SM-like couplings. Conversely, in the case of production (right plot) the EW corrections reach a comparable value only in the strip of very low , and they quickly drop below as soon as  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 GeV [20], while in the case of and 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 and . 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 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 , , 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 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  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 production in the decoupling limit never goes beyond . For what concerns production, the effect of the interference between quark and squark contributions becomes significant only for very large and moderate , where the gluon-fusion cross section is tiny anyway. The largest effect, a suppression by , is found in the light-stau scenario for  GeV and , where the cross section is of the order of a tenth of a femtobarn. The SUSY contributions to production, already small in the light-stop scenario because they only arise at two loops, are further suppressed in the , , 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 and , the radiative corrections to the element of the CP-even Higgs mass matrix suppress significantly the mixing angle , so that the coupling of to taus – which is proportional to – 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 to bottom quarks is suppressed. As a result, in the tau-phobic scenario the behavior in the 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 production shows some enhancement with even for large values of , while for small the total cross section for production has a milder dependence on than in the other scenarios. Also, the suppression of the coupling to bottom quarks makes the contribution of bottom-quark annihilation to production smaller than in the other scenarios. Finally, the tau-phobic scenario is characterized by third-generation squark masses around  TeV, and by a value of the superpotential Higgs-mass parameter,  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 , , and light-stau scenarios, and in the case of pseudoscalar production it is even larger than in the light-stop scenario.

## 3 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 , 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, -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.

### 3.1 Scale dependence of the cross section

In this section we study the dependence of the cross section for Higgs production on the renormalization scale at which the relevant couplings in the partonic cross section are expressed, and on the factorization scale 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 and , 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 and 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, 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 . In our study of gluon fusion we take 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 following the LHC-HXSWG prescription [3]: we consider seven combinations of renormalization and factorization scales, defined as the set of the pairs obtainable from the two sets and , with the additional constraint that (i.e., we treat the variations of the ratio 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 ,

 σ− ≡ min(μR,μF)∈Cμ{σ(μR,μF)} ,       σ+ ≡ max(μR,μF)∈Cμ{σ(μR,μF)} , (2)

and define the relative scale uncertainty of the cross section as , where

 Δ+μ ≡ σ+−σ(¯μR,¯μF)σ(¯μR,¯μF)  ,       Δ−μ ≡ σ−−σ(¯μR,¯μF)σ(¯μR,¯μF)  . (3)

In figures 8 and 9 we show the contours of equal for scalar and pseudoscalar production in the 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 (i.e., NLO) and in the VHML at (i.e, NNLO). Its residual scale dependence amounts to an effect, with the exception of some mass-dependent effects at , which are known to be numerically small [14]. The bottom and sbottom contributions are included at the NLO and they account for an effect. The stop contributions are included through the NNLO, see section 3.4, but their effect on scale dependence is also of because we neglect the genuine three-loop terms. Finally, while the EW corrections are computed at , their inclusion as a fully factorized term at the NLO causes their effect on scale variation to be of , 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 plane. The uncertainty is lowest, around , where the top contribution dominates: this is the case for production (left plot in figure 8) in the decoupling region, where the uncertainty stabilizes to roughly at large (i.e., slightly smaller than the we obtain for the same Higgs mass in the SM); for production (right plot in figure 8) in the strip with  GeV, as well as when and  GeV; for production (figure 9) in the strip with . In contrast, the scale uncertainty exceeds in the regions where the bottom contribution is enhanced or downright dominant: at large for and production, and at small for production.

The plots for and production in figures 8 and 9 show additional structures. In the case of production, the scale uncertainty becomes very large for and  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 production, the structure visible for  GeV is associated to the cusp-like behavior of the top contribution to the gluon-fusion amplitude around the threshold . Another feature of and 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 production, the scale dependence in the decoupling region is similar to, or even bigger than, the one in the SM. For 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 and  GeV occur at higher values of .

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 and , 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 and , with the additional constraint that (again, we discard the two pairs with the largest variation of around its central value, which in this case is ). 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 . This is due to the fact that, in eq. (3), the -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 and .

In figures 10 and 11 we show the scale dependence of the cross section for scalar and pseudoscalar production, respectively, as a function of in the light-stop scenario with . In the upper part of each plot, the solid line denotes the cross section for bottom-quark annihilation computed with the central scale choice , 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 production can be as large as for low values of , then it stabilizes to roughly in the decoupling region where becomes independent of . In contrast, the relative scale uncertainty of the cross section for the production of (right plot in figure 10) and (figure 11) decreases as (and hence ) 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.

### 3.2 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, (here GeV). In the MSSM the tree-level couplings of the neutral Higgs bosons to bottom quarks are modified as follows:

 Yhb = −sinαcosβ YSMb ,       YHb = cosαcosβ YSMb ,       YAb = tanβ YSMb , (4)

where is the mixing angle in the CP-even Higgs sector. In the decoupling limit, , the mixing angle simplifies to , so that the coupling of to bottom quarks is SM-like, while the couplings of and are both enhanced by .

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 -enhanced SUSY contributions are resummed in effective Higgs-bottom couplings.

#### Scheme and scale dependence of the bottom mass

The parameter enters the expression for the gluon-fusion amplitude with two distinct roles: as the actual mass of the bottom quarks running in the loops, and as a proxy for the Higgs-bottom coupling , where . The numerical value of depends strongly on the renormalization scheme and scale: an mass  GeV corresponds to a pole mass  GeV at three-loop level, whereas evolving up to a scale of the order of the typical energy of the gluon-fusion process decreases significantly its value. For example, if we evolve at four-loop level the bottom mass up to the scale at which we express the strong gauge coupling, , we obtain  GeV for  GeV. While any change in the definition of the bottom mass and Yukawa coupling entering the one-loop part of the amplitude is formally compensated for, up to higher orders, by counterterm contributions in the two-loop part, the numerical impact of such strong variations on the prediction for the gluon-fusion cross section can be significant.

To illustrate this point, we identify the mass of the bottom quarks in the loops with the pole mass , and consider the dependence of the gluon-fusion cross section on the prescription for the Higgs-bottom coupling , focusing on . In the light-stop scenario with  GeV and , where both Higgs scalars are relatively light and have enhanced couplings to the bottom quark, the effect of extracting from the mass instead of the pole mass leads to a decrease in the cross section for production, and a decrease in the cross section for production. The use of would instead decrease the cross section for production by , and the one for production by , with respect to the values obtained with . As a second example, we take the light-stop scenario with  GeV and , where the lightest scalar has SM-like couplings to quarks. In this case the cross section for production varies by less than when choosing among the three options discussed above for the definition of . For the heaviest scalar , on the other hand, the changes in the cross section relative to the value derived with amount to and when is extracted from and , 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 to or to might seem preferable to the one of using , 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 and the bottom Yukawa coupling with , where is a generic renormalization scale, the contribution of diagrams with bottom quarks and gluons to the two-loop part of the amplitude reads

 A2ℓb(τ)  ∝  CF[FCF(τ)+F1ℓ1/2(τ)(1−34lnm2bμ2b)]+ CAFCA(τ) , (5)

where and are color factors, , and we omit an overall multiplicative factor. Truncating the functions at the first order in an expansion in powers of , one finds [11]

 F1ℓ1/2(τ) = −2τ(1−14L2bϕ)+ O(τ2) , (6) FCF(τ) = −τ[5+95ζ22−ζ3−(3+ζ2+4ζ3)Lbϕ+ζ2L2bϕ+14L3bϕ+148L4bϕ]+ O(τ2) , (7) FCA(τ) = −τ[3−85ζ22−3ζ3+3ζ3Lbϕ−14(1+2ζ2)L2bϕ−148L4bϕ]+ O(τ2) , (8)

with

 Lbϕ ≡ ln(−4/τ) = ln(m2ϕ/m2b)−iπ . (9)

The equations above show that the two-loop bottom contribution to the gluon-fusion amplitude contains powers of , and that the choice does eliminate some of the logarithmically enhanced terms. Similarly, relating the coupling entering the one-loop part of the amplitude to the pole mass eliminates the whole piece proportional to in eq. (5). Each of the two remaining terms, and , also contains powers of , but for realistic values of 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 , which originates from diagrams with three- and four-gluon interactions. In that case no significant cancellation occurs, and the amplitude is not minimized when is extracted from or . 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 .

In the case of the Higgs coupling to photons, the problems related to the ambiguity in the definition of have been solved with a resummation of the leading and next-to-leading logarithms of the ratio  [89]. 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 , and we consider the difference between the results obtained using and those obtained using 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  GeV, this procedure – also advocated by the LHC-HXSWG in ref. [3] – results in an uncertainty of 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 in the regions of the 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 when is related to . 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 -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

 ˜Yhb = Yhb1+Δb(1−Δbcotαtanβ),   ˜YHb = YHb1+Δb(1+Δbtanαtanβ),   ˜YAb = YAb1+Δb(1−Δbcot2β), (10)

where are the tree-level Higgs-bottom couplings defined in eq. (4), and, retaining only the contribution from diagrams with sbottoms and gluinos, the -enhanced term reads

 Δb = 2αs3π m~gμtanβm2~b1−m2~b2⎛⎝m2~b1m2~b1−m2~glnm2~b1m2~g−m2~b2m2~b2−m2~glnm2~b2m2~g⎞⎠. (11)

In the limit , where , the superparticle contributions encoded in 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 .

In refs. [39, 40] it was shown that, in the calculation of processes that involve the Higgs-bottom couplings, the -enhanced corrections can be resummed to all orders in the expansion in powers of 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 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 in the resummed one-loop amplitude with the term of the expansion of in powers of . Depending on the choice of renormalization scheme for the parameters in the sbottom sector, additional -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 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 or in the two-loop part is formally of higher order, i.e., it amounts to three-loop terms that are suppressed by a factor with respect to the dominant three-loop terms of 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 . We found that this choice allows us to reproduce – after an expansion in powers of – the three-loop result that can be inferred from ref. [40], where the sub-dominant terms proportional to were also resummed in the effective couplings.

For large values of , the factor 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 . For positive the correction suppresses the couplings, reducing the overall relevance of the bottom contribution to gluon fusion. On the other hand, for negative the correction enhances the couplings, which diverge as approaches . As a consequence, when is large and negative the result for the gluon-fusion cross section is extremely sensitive to the precise value of , 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 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 contribution in eq. (11). In our numerical analysis we use by default the full one-loop result for as computed by FeynHiggs, which allows us to resum in our prediction for the Higgs-production cross section also the -enhanced corrections of electroweak origin.

Another improvement in the calculation would come from the inclusion of the dominant two-loop contributions to , 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 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 vary by roughly 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 by varying by the result provided by FeynHiggs.

In general, the impact of the uncertainty of 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 discussed in the previous section, the 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  GeV and , where both Higgs scalars have enhanced couplings to bottom quarks. The superpotential parameter has positive sign, and the corrections suppress the effective couplings . We find that the cross sections for and production in gluon fusion increase by and , respectively, if the value of is reduced by , while they decrease by and , respectively, if is increased by . The effect is larger if is taken negative, so that the corrections enhance the effective couplings. In that case the dependence on is reversed: if we consider the same point in the light-stop scenario but flip the sign of , the cross sections for and production in gluon fusion decrease by and , respectively, when is reduced by , while they increase by and , respectively, when is increased by .

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 . Also in this case, we can estimate the uncertainty by varying by the value of provided by FeynHiggs.

### 3.3 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 bottom-mass 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 different PDF sets (replicas), all equivalent from the statistical point of view in the description of the data. Any observable has to be computed 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 at the LO and is subject to very large QCD corrections, of , at the NLO. Each PDF collaboration recommends a different central value for , generating a spread of the central predictions for the Higgs-production cross section. The combination of the PDF and uncertainties (henceforth, PDF+) 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+ uncertainty bands of the three collaborations, where for each group the preferred central value is adopted [94]. Following this reference we take 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+ 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+ 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+ uncertainty is at the sub-percent level.

To assess the PDF+ 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+ 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.9 We find that the uncertainty has an almost constant behavior when the mass of the produced Higgs boson is lighter than  GeV, and that it increases for larger mass values: for example, at the NNLO, the PDF+ uncertainty of the cross section for bottom-quark annihilation amounts to for  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  GeV. Converted to the mass via a three-loop QCD calculation, this corresponds to  GeV, which differs both from the value recommended by the LHC-HXSWG,  GeV [57, 59], and from the current PDG value,  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 ranging from