Hadronic Higgs production through NLO + PS in the SM , the 2 HDM and the MSSM

The next-to-leading order (NLO) cross section of the gluon fusion process is matched to parton showers in the MC@NLO approach. We work in the framework of MadGraph5_aMC@NLO and document the inclusion of the full quark-mass dependence in the Standard Model (SM) as well as the state-of-the-art squark and gluino effects within the Minimal Supersymmetric SM embodied in the program SusHi. The combination of the two programs is realized by a script which is publicly available and whose usage is detailed. We discuss the input cards and the relevant parameter switches. One of our focuses is on the shower scale which is specifically important for gluon-induced Higgs production, particularly in models with enhanced Higgs-bottom Yukawa coupling.


Introduction
Higgs production proceeds predominantly through gluon fusion in a large number of theories, including the Standard Model (SM). The recently discovered resonance [1,2] in searches for a Higgs boson is fully consistent with the SM picture, 1 so far. Still, the measured Higgs boson may be embedded in an enlarged Higgs sector with respect to the one of the SM which predicts only a single physical particle breaking the electro-weak symmetry. Two-Higgs doublet models (2HDM's) such as the Minimal Supersymmetric SM (MSSM) are among the most popular theories with enlarged Higgs sectors. Such theories inevitably require the existence of further physical Higgs particles. A 2HDM predicts three neutral Higgs bosons: a light (h) and a heavy 1 See Refs. [3][4][5] for a theoretical overview. a e-mail: hmantler@cern.ch b e-mail: mariusw@physik.uzh.ch (H ) scalar, and a pseudo-scalar (A); as well as two charged Higgs particles (H ± ). Almost all 2HDM and MSSM scenarios that are in agreement with the experimental bounds feature a light scalar which is SM-like in its couplings to vector bosons and fermions, while the other Higgs bosons are heavier and, therefore, escaped detection up to now. Indeed, the experimental search for other Higgs resonances is one of the major focuses regarding the discovery of physics beyond the SM (BSM) in the second run of the Large Hadron Collider (LHC).
Higgs production through gluon fusion is mediated by a colored particle. In the SM, the top quark gives the dominant contribution to the cross section [6][7][8]. While also the bottom quark gives a sizable contribution, the effects due to other quarks are small and therefore usually neglected. In the 2HDM and the MSSM the Higgs-bottom Yukawa coupling can be enhanced with respect to the one of the top quark and the bottom loop may even constitute the dominant contribution to the cross section. In those models it is stringently required to include the bottom-quark contribution. 2 The gluon fusion cross section is known at the next-toleading order (NLO) in the SM including top-and bottommass effects [18,19], in the 2HDM and in the MSSM including contributions from squarks and gluinos [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. For the top quark, an effective field theory approach can be applied in which the top quark is considered to be infinitely heavy and can be integrated out from the full theory. In this approximation, Higgs production has been calculated up to the next-to-NLO (NNLO) inclusively [35][36][37] as well as fully differential [38][39][40]. Electro-weak contributions and effects beyond NNLO in the heavy-top approximation have been studied in Refs. [41][42][43][44][45][46][47][48][49][50] for example, while there was a large effort [51,52] to push the accuracy to next-to-NNLO ( N 3 LO) which has been succeeded very recently [53]. Finite top-mass effects have been shown to be small for both the inclusive cross section at the NNLO [54][55][56][57][58][59] and differential quantities [60,61] as long as no kinematical scale (such as the transverse momentum of a particle) that is not integrated out exceeds the top-mass threshold.
The full dependence of the top-and the bottom-mass at the NLO has been included so far in a POWHEGtype [62] matching to parton showers (PSs) [33], the analytically resummed transverse momentum spectrum of the Higgs boson at NLO + NLL [63], a MC@NLOtype [64] matching to the Herwig Monte Carlos [65][66][67], the NNLO + NNLL jet-vetoed [68] and the fully differential NNLO [69] cross section; and in some approximated form recently also in the NNLOPS approach [70,71]. Furthermore, the 2HDM as well as supersymmetric effects from squarks and gluinos within the MSSM [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] have been implemented in the first two approaches from that list [33,72]. In this manuscript, we report on a new implementation of NLO QCD corrections in the SM, 2HDM and MSSM applying the MC@NLO-type matching to both Herwig and Pythia showers. We work in the framework of MadGraph5_aMC@NLO [73] and combine its capabilities with the corresponding amplitudes provided by SusHi [74]. The linking of SusHi to MadGraph5_aMC@NLO is realized by a script 3 aMCSusHi. Its usage as well as the application of the combined code to obtain cross section predictions in the SM, the 2HDM and the MSSM is detailed in this paper. 4 The manuscript is organized as follows: In Sect. 2 we present a brief overview of the elements of the computation at hand. Section 3 is dedicated to introduce aMCSusHi and is separated in three parts which cover: how to use the script (Sect. 3.1), how to run the resulting code (Sect. 3.2) and how to treat the shower scale (Sect. 3.3). We will show a brief application of the code to phenomenological results in Sect. 4 and conclude in Sect. 5.

Outline of the calculation
The goal of this paper is to present a tool which allows for the computation of arbitrary infra-red safe differential observables at both the parton and the hadron level for the production of neutral Higgs bosons via gluon fusion in the SM, the 2HDM and the MSSM by matching the NLO cross section to a shower.
The relevant NLO matrix elements are taken from Ref. [74], which include both SM-like contributions and sbottom, stop and gluino effects. Examples of corresponding Feynman diagrams are illustrated in Fig. 1. They are combined and matched to a parton shower by the well-known MC@NLOmethod. The matched cross section in the MC@NLO framework can be written symbolically as  Fig. 1a-c. The NLO virtual and real corrections are governed by diagrams like the ones shown in Fig. 1d-g and h, i, respectively, and similar ones with quark loops replaced by squark loops.
Equation (1) is implemented for all standard parton showers [76][77][78][79][80][81] in the fully automated framework MadGraph5_aMC@NLO. This code determines NLO QCD corrections to arbitrary scattering processes at the LHC. On the basis of UFO models [82], the code even allows one to carry out computations in any theory beyond the SM in a general manner as soon as the renormalization is known and implemented in a UFO model; see Refs. [73,[83][84][85] for further information. However, the Higgs production mode through gluon fusion is special, being loopinduced already at the LO. Such processes cannot be handled in a fully automated manner by any code to date, since it requires the automation of two-loop amplitudes which is beyond current technology. Therefore, we have treated Higgs production through gluon fusion in the SM, 2HDM and MSSM as a special case, by linking the relevant amplitudes from SusHi. Furthermore, as far as the MSSM is concerned SusHi requires a link to FeynHiggs [86][87][88][89][90][91][92][93][94][95][96] which evaluates the corresponding Higgs masses and couplings in user defined scenarios. Setting up the SusHi amplitudes in MadGraph5_aMC@NLO is handled by a publicly available script called aMCSusHi, which is automated to create the gg → φ process folder; download SusHi and FeynHiggs; compile, install and link them; and replace the relevant amplitudes in the process folder. In the upcoming section, we describe the application of the script and explain the necessary steps to obtain phenomenological results. Fig. 1 A sample of Feynman diagrams for gg → φ contributing to the NLO cross section; a-c LO, d-g virtual, and h, i real corrections. The graphical notation for the lines is: solid straight = quark; spiraled = gluon; dashed = scalar (squark or Higgs); This code is based on MadGraph5_aMC@NLO and SusHi. For further information on theses codes we refer the reader to corresponding publications [73,74]. After using the script to set up the code, we will focus on the relevant user inputs to obtain phenomenological predictions for Higgs cross sections in the SM, the 2HDM and the MSSM. The first argument is mandatory and determines the path to the process folder for gg → φ that is generated by the script. The only requirement is that this folder has to be defined as a sub-folder of the MadGraph5_aMC@NLO directory. The second and third arguments are optional. If there are compiled versions of FeynHiggs and SusHi available on your computer you can give the names of the folders that contain the files libFH.a and libsushi.a, respectively. When executed with two (one) argument(s) the script will ask whether it should automatically download and install SusHi (and FeynHiggs). The script will always download the latest versions of these codes. While running, the script requires some user inputs: It asks whether or not SusHi (and FeynHiggs) should be downloaded, in which folder they should be installed (default is inside the ggH-folder ) or where to find SusHi (and FeynHiggs) if already installed. The user is simply required to follow these on-screen instructions. Furthermore, the script creates log-files in the working directory for the download ("XX_curl.log"), the configure command ("XX_conf.log") and the compilation ("XX_make.log"), where XX = FH for FeynHiggs and XX = SusHi for SusHi. These files are supposed to give complementary information for any kind of troubleshoot-ing. 6 Further information about the aMCSusHi script provides the README file.

Running the code
Once the gg → φ process folder has been set up by the script, the run can be started directly from the ggH-folder by typing > ./bin/generate events and following the usual steps in MadGraph5_aMC@NLO to choose the run-modes (order, shower or fixed-order, madspin). Before running the code, one may want to modify the input settings. In the following we will discuss the differences between aMCSusHi and the ordinary MadGraph5_aMC@NLO code regarding the input files.
The input cards can be found under ggH-folder /Cards/, where the param_card.dat, run_card.dat and shower_card.dat contain all the essential information. The run card controls the usual parameters, e.g. the renormalization (μ R ) and factorization scale (μ F ). Note that by default the flags fixed_ren_scale and fixed_fac_scale are set to false, so that these two scales are chosen on an event-wise basis. Their values are specified in the last routine of ggHfolder /SubProcesses/setscales.f which, due to the default option dynamical_scale_choice = −1 in the run card, sets where i runs over all final state particles and m i and p T (i) are their mass and transverse momentum, respectively. This choice is reasonable, since it respects effects from hard radiation and corresponds to a value of m φ /2 in the soft/collinear limit which is the current recommendation for the total inclusive gg → φ cross section [3,5].
Also the shower card in aMCSusHi contains no new information and has the usual functionality. Since MadGraph5_aMC@NLO supports all standard parton showers, for the first time Pythia6 and Pythia8 can be applied at NLO + PS to SM Higgs production in the full theory in the MC@NLO framework. In general, it is advisable to apply the most recent versions of the showers for meaningful physics runs.
The param_card.dat, on the other hand, receives some significant changes by the aMCSusHi script. The new version basically combines the orignial parameter card from MadGraph5_aMC@NLO with the input file from SusHi which are both written in the SUSY Les Houches accord (SLHA) format [97] and, therefore, easily connectable. We will address the different options in the param_card.dat in more detail, since there are a number of changes and some of the original parameters lose their functionality. In the SLHA format inputs are organized in blocks which have different entries that are characterized by a number. For simplicity, we introduce the following shorthand notation: Block example[i] corresponds to entry i in Block example. E.g., entry 25 of Block mass (Block mass [25]) in the SLHA format is devoted to the Higgs mass in the SM, which is required as an input in the param_card.dat. A typical parameter card of aMCSusHi in the SM is shown below: Most of the inputs are self-explanatory due to the comments initiated by the hash symbol # after the entries. Furthermore, the standard SLHA blocks match the universal convention of Ref. [97]. Some of the inputs, though, require further comments. In the Block mass all parameters have the expected function, except for the top and the bottom mass, Block mass [6] and Block mass [5], respectively. While the former only affects and is required for the shower, the latter can be omitted completely. Instead, due to the link to SusHi the top mass that is used for the top loop and the top Yukawa is set in Block sminputs [6] and is expected to be on-shell. For the bottom mass, on the other hand, SusHi allows for three different choices: on-shell scheme or MS scheme with m b (m b ) or m b (μ R ), which can be switched in Block renormbot [1] by a value between 1 and 3. 7 Also here the recommendation is to use the on-shell scheme which according to Refs. [18,98] ensures the cancelation of large logarithms ln(m h /m b ) at NLO QCD, while the MS scheme does not, due to an incomplete resummation of these terms. The on-shell m b value is determined by Block renormbot [4], while when the MS scheme is chosen the corresponding input of m b (m b ) is set in Block sminputs [5]. The other entries of Block sminputs again have the same impact as in the usual MadGraph5_aMC@NLO code. The same is true for Block yukawa. For the decay of light and heavy Higgs bosons one may specify a finite width of the Higgs boson in the respective BSM scenario by using Decay 25 irrespective of whether the light or the heavy Higgs boson is considered. At this point we shall remark that the particle identification number in the generated event files is always 25 regardless of the Higgs boson under consideration. 8 This is irrelevant for the production (which is correctly computed through the SusHi amplitudes), but it plays a role for the decay where the shower will consider particle 25 to be the light Higgs, which is indeed fine for any scalar Higgs, but a problem for pseudo-scalar ones. Therefore, decays of a pseudo-scalar Higgs are currently not supported in the official version of aMCSusHi. A user interested in decaying the pseudo-scalar Higgs is strongly encouraged to contact us.
The other parameters are relevant to SusHi. Block sushi [1] 8 Bear in mind that this has to be taken into account to identify the Higgs in the analysis of the showered events. 9 A link to 2HDMC [99] with the corresponding input convention is currently not supported.
Additionally to the inputs which we defined already for the SM the following parameters have to be set in the 2HDM: Block renormbot [2] specifies whether or not a resummation 10 of terms enhanced by tan(β) is applied through reweighting of the bottom Yukawa coupling as described in Ref. [74]; Block 2hdm determines which type of the 2HDM is used; the value of tan(β) is set through Block minpar [3]; and the mixing angle α corresponds to the entry in Block alpha.
The computation of MSSM Higgs cross sections requires Block extpar, Block feynhiggs and Block renormsbot in addition, which fix the parameters of the third family of quarks and squarks, determine the FeynHiggs inputs and yield information on the renormalization of the sbottom section, respectively. We will not provide any further information on these blocks, instead, we refer to the SusHi manual [106] and the FeynHiggs man pages [107]. Moreover, Block alpha can be omitted in the MSSM and the Higgs masses in Block mass have no effect, since they are determined by FeynHiggs, once Block feynhiggs is present. The MSSM Higgs mass that has been computed and applied in the run is provided to the user in Block mass [25] of the parameter card, which will be overwritten by the mass of the respective Higgs boson at the beginning of each MSSM run.
So far we did not comment on the Block factors. It allows one to turn on and off individual contributions in all models. In fact, it even provides the possibility to rescale the respective Yukawa couplings by choosing values different from 0 and 1. With Block factors [1] one can include the charm quark in the computation. This requires one to specify its MS mass m c (m c ) in Block sminputs [8] which is then translated to its on-shell mass. Furthermore, Block factors [2] and Block factors [3] multiply the top and the bottom Yukawa, respectively. In the MSSM, the stop Yukawa is rescaled by Block factors [4] and the sbottom one by Block factors [5].
For further information on the input cards we refer the reader to Ref. [73] of MadGraph5_aMC@NLO, the manual of SusHi [106] and the man pages of FeynHiggs [107]. Three example parameter cards are provided in the folder ggH-folder /Cards; one for the SM (param_card.dat_SM), the 2HDM (param_card.dat_2HDM_scenB) and the MSSM (param_card.dat_MSSM_mhmodp). They match the scenarios that we study in the result section of this paper.

Choosing different shower scales
The choice of the shower scale is a very peculiar one in the gluon fusion process. In presence of the bottom-quark loop, factorization of soft and collinear radiation maybe spoiled at scales significantly smaller than the Higgs boson mass. This was pointed out by Ref. [69] in the context of analytic transverse momentum resummation. On the other hand, these terms might well be treated as a finite remainder as long as their impact remains moderate [68].
Due to their additive matching of the resummed lowp T region with the fixed-order distribution valid at large transverse momenta, analytic p T -resummation and the MC@NLO-method are quite similar. In both cases there is a scale associated with that matching, the resummation scale Q res and the shower scale Q sh , respectively. These scales can be interpreted as transition scales that separate the soft/collinear from the hard physics, very similar to the factorization scale of the PDFs. In other words, they define the range where resummation, and therefore the shower in MC@NLO, takes effect. Their value has to be chosen of the order of the typical scale of the problem.
In gluon fusion, the typical scale depends on the quark considered in the loop. Since m t ∼ m φ , there exist only two relevant scales for the top-quark loop (m φ and p T ) and the shower scale can be chosen of the order of the Higgs mass. When considering the bottom loop, on the other hand, we face a three-scale problem (m φ , m b , and p T ) which has not been solved to date. However, it has been suggested [69] to apply a lower transition scale to the bottom contribution, which, in particular, respects the fact that soft/collinear factorization is valid only up to smaller scales for the bottom loop. In Ref. [72] it was further proposed to separate three contributions according to their Yukawa couplings: the square of the top and the bottom, and their interference; and choose separate shower/resummation scales for all of them. This splitting allows for a model independent treatment of the problem by a rescaling of the individual contributions with the respective top and bottom Yukawas of a specific scenario in the 2HDM as well as the MSSM when neglecting squark effects. In the literature, two pragmatic approaches with physical motivation have been presented [72,108] to determine separate scales for the three contributions. Their comparison will be studied elsewhere [109]. When studying phenomenological results in Sect. 4 we will apply the scales from Ref. [72] (referred to as "HMW" in what follows).
The separation of the bottom contribution (including the interference) from the top one with different shower scales (Q b and Q t , respectively) requires three runs in MadGraph5_aMC@NLO, which have to be combined as follows: To obtain all three contributions of different Yukawa origin with different scales that allows for a model independent treatment, on the other hand, MadGraph5_aMC@NLO has to be run five times: The scales Q t , Q b , and Q tb determine the scale for top, the bottom and their interference, respectively. As indicated before, the individual contributions can be separated using the Block factors in the parameter card. The shower scale in MadGraph5_aMC@NLO cannot simply be accessed through the input cards, since it requires an advanced user to be familiar with its specific treatment in the code. MadGraph5_aMC@NLO does not use a simple fixed scale for Q sh , instead, it statistically extracts the shower scale from a distribution peaked at a specific value. The user can only change the range of the interval of the distribution which of course also affects the peak. Therefore, we identify Q t , Q b , and Q tb in Eqs. (2) and (3) with the peak of the respective shower scale distributions.
The so-called shape parameters define the interval of the distribution from which the shower scale is picked on an event-wise basis. They can be specified in the include file ggH-folder /SubProcesses/madfks_mcatnlo.inc, where the relevant part is given by (MadGraph5_aMC@NLO default values): The parameters frac_low and frac_upp are used to compute the upper and lower bounds of the Q sh distribution, which will be explained in more detail below, while scaleMClow allows one to set an absolute value of the lower bound on Q sh and scaleMCdelta is used to apply a minimal value to the size of the distribution interval. In formulas the interval is defined by 11 where s 0 is the Born-level partonic center of mass energy squared. Evidently, scaleMClow and scaleMCdelta only take effect if the interval obtained through frac_low and frac_upp does not meet the corresponding restrictions. The corresponding Q sh distribution is peaked around For a 2 → 1 process like gluon fusion this relation is an identity and √ s 0 equals the mass of the final state particle, i. e., the Higgs mass in the case of gluon fusion. To change the peak-value to its half, e. g., for shower scale variations, one can simply divide frac_low and frac_upp by a factor of two. In this sense, it is convenient to keep the ratio between frac_low and frac_upp a constant, which in the default setup of MadGraph5_aMC@NLO is a factor of ten. Under this prerequisite, in order to choose a specific shower scale Q sh for gg → φ we simply have to determine Here and in what follows, we associate Q peak with the shower scale Q sh and vice versa unless stated otherwise. After modifying the corresponding include file accordingly, 11 For further information we refer the reader to Sect. 2.4.4 of Ref. [73].
MadGraph5_aMC@NLO has to be recompiled. This can be achieved by typing > make clean inside the ggH-folder , which forces a recompilation during the next run of the code.
In Sect. 4 we show some applications of the aMCSusHi code and study the effect of different treatments of the shower scales.

Results: brief application
The gg → φ process folder created by the aMCSusHi script preserves all the highly convenient features that come with MadGraph5_aMC@NLO. Besides many others, this entails an interface to the most common showers, the fully automatic determination of scale and PDF variations without any extracosts of computing time [110], the creation of any number of observables with a single run and analysis routines available for the most important processes including the gluon fusion Higgs production mode.
aMCSusHi allows one to compute gluon-induced Higgs production including the complete dependence on the quark masses in the SM for the first time in a MC@NLO-type matching applying all versions of Pythia and Herwig Monte Carlos. While previous computations did only feature the Herwig showers [65][66][67], phenomenological results exist to our knowledge only for Herwig6 [111]. As a first application we therefore study the impact of different showers on the top-and bottom-mass effects with respect to the heavy-top approximation at the 13 TeV LHC. For this purpose, Fig. 2 shows the ratio of the NLO + PS computation including mass effects and the corresponding cross section in the heavy-top limit as a function of the transverse momentum of the Higgs for different Monte Carlos: Pythia8 (black, solid), Herwig++ (red, dotted), Pythia6 p T -ordered (blue, dashed with points), Pythia6 Q-ordered (green, dash-dotted with open boxes) and Herwig6 (yellow, solid with filled boxes). We apply the MSTW2008 68 % CL NLO PDF set [112] with the corresponding value of the strong coupling constant. The shower scale has been chosen as Q sh = m h /2 in all cases, while for μ F and μ R we use the defaults specified in Sect. 3.2. Clearly, the mass effects are hardly dependent on the specific Monte Carlo, which is particularly evident at small ( p T 50 GeV) and large ( p T 150 GeV) transverse momenta. Nevertheless, there are some visible differences in the intermediate region which consistently discriminate the Herwig from the Pythia showers. Overall, they are at most 5 % though and therefore still moderate. Figure 3 shows the effects of quark masses with respect to the heavy-top approximation as well, but for different choices of the associated shower scales. In all cases, the denominator and therefore the distribution in the heavy-top limit is computed with the respective scale of the top contribution Q sh = Q t . As we observed before the Monte Carlo dependence is quite small; therefore, we only consider the Pythia8 shower. For reference the black solid curve is the same as in Fig. 2 with Q sh = m h /2 for all contributions. We compare it to the scales choices proposed in Ref. [69], which imply setting the shower scale of the bottom contribution (including the interference) to the bottom mass following Eq. (2) (red dotted curve). For the blue dashed curve with points we chose the HMW scales determined in Ref. [72] which can be found in Table 1 of that paper, applying a three-scale approach according to Eq. (3) (Q t = 49 GeV, Q tb = 34 GeV and Q b = 23 GeV). The green dash-dotted curve with open boxes serves mostly for comparison with previous Herwig6 results [111] which were computed with the scales of Ref. [69] as well.
For the p T distribution in Fig. 3a, the change of the scale of the bottom contribution to Q b = m b has a significant impact on the mass effects at small and intermediate transverse momenta. It develops an extremely steep drop at small transverse momenta which due to unitarity affects also the intermediate p T -range in the opposite direction. The benefit of the usage of such a low scale is clearly disputable. While the Herwig6 curve agrees rather well with previous result of Ref. [111] becoming flat for p T 5 GeV, the Pythia8 curve develops a steep increase in this region. This signals a significant Monte Carlo dependence at very small p T which is not observed for larger Q b scales. Furthermore, the rigorously low value also poses a technical problem in the code regarding the fact that the default shower scale choice in MadGraph5_aMC@NLO, as explained in Sect. 3.3, is a distribution. In order to solve this problem, we had to use a fixed value of Q sh = Q b by setting frac_low = frac_upp = Q b /m h and scaleMClow = scaleMCdelta = 0.
Considering the HMW scales, the scale of the bottom contribution is not chosen at such low values. We find that the mass effects in this case (blue dashed line with points) are rather similar to the ones where all scales are set to m h /2 (black solid line), although the individual HMW scales being quite different from this value. Looking at the rapidity distribution in Fig. 3b, on the other hand, we observe the expected feature of being essentially insensitive to any choice of the respective shower scales. We shall note at this point that simply due to their inclusion in the default analysis we were able to produce a large number of further observables at no additional computing cost.
To demonstrate the range of applicability of aMCSusHi, we consider two realistic BSM scenarios in Figs. 4 and 5: the heavy Higgs boson in Scenario B of Ref. [113] (a bottom dominated 2HDM scenario); and the pseudo-scalar Higgs boson in the m mod+ h (800, 40) MSSM scenario [114] defined in Table 2 of Ref. [72]. The corresponding input files can be found in the folder ggH-folder /Cards.
In Fig. 4 we study the transverse momentum distributions of the Higgs and the hardest jet, while Fig. 5 depicts their rapidity distributions. In both cases we apply the HMW scales of Ref. [72]. At low transverse momenta the p T distributions have similar shapes comparing the red dotted (LO + PS) to the black solid curves (NLO + PS). This can easily be inferred from the first inset where all curves are normalized to the black solid line in the main frame. However, it is well known that at the LO + PS p T distributions yield unphysical results for transverse momenta beyond the shower scales indicated by a steep drop. Note that both curves are normalized to the same (the NLO) cross section. Continuing the comparison at hand, we observe a significant reduction of the scale uncertainties shown in the lower inset, where the bands are obtained by dividing the upper and  4 Transverse momentum distribution of a the heavy Higgs boson and b the associated hardest jet computed in a bottom dominated scenario of the 2HDM (see text for details). The graphical notation is the following: the black solid curve shows Pythia8 at NLO + PS, the red dotted curve is the same at LO + PS (normalized the NLO) and the blue dashed one with points corresponds to the fixed-order curve at NLO lower bound of the respective cross section by the same central cross section as in the first inset. The uncertainties correspond to the independent variation of all unphysical scales (μ R , μ F , Q t , Q b , Q tb ) by a factor of two. Comparing NLO + PS to the NLO fixed-order result denoted by fNLO, we observe the expected matching toward large transverse momenta.
In order to compare shapes, the rapidity distributions in Fig. 5 are normalized in a way that the sum of their bins yields 1. We see that for the Higgs rapidity in Fig. 5a all curves agree extremely well in terms of shape up to the forward region in which, nevertheless, the deviations are still well within the respective uncertainty bands. 12 For the rapidity distribution of the hardest jet the same is true for the two showered results, while the fNLO distribution, on the other hand, agrees only scenario [114] of the MSSM with M A = 800 GeV und tan β = 40. The graphical notation is the same as in Fig. 4. All curves are normalized so that their bins add up to 1 in the central region |y( j 1 )| 3, but it features a significant enhancement when the hardest jet is more forward. In this region the cross section will receive large effects of collinear radiation which renders the shower to yield the more reliable description.

Conclusions
In this article we presented the new tool aMCSusHi which is a link between MadGraph5_aMC@NLO and SusHi for the computation of Higgs cross sections in gluon fusion at hadron colliders. The code gives NLO + PS accurate results in the SM, 2HDM and MSSM. The inputs in the MSSM are conveniently controlled through a link to FeynHiggs. We discussed the specific treatment of the shower scale in MadGraph5_aMC@NLO and pointed out its special role in the context of gluon-induced Higgs production. In the phenomenological part we study the impact of different shower scale choices on the mass effects in the SM. Furthermore, we studied results for 2HDM and MSSM Higgs production as an application of aMCSusHi. The aMCSusHi script can be downloaded from https://cp3.irmp.ucl.ac.be/projects/ madgraph/wiki/aMCSushi.
As an outlook, one may improve the SM prediction for the Higgs production mode through gluon fusion by merging the NLO + PS cross section for gg → h in the full theory with higher multiplicities computed in the heavy-top effec-tive field theory. This is certainly beyond the scope of the present paper and is left for a future publication.