Heavy charged Higgs boson production at the LHC

In this paper we study the production of a heavy charged Higgs boson in association with heavy quarks at the LHC in the two-Higgs-doublet model. We present for the first time fully-differential results obtained in the four-flavour scheme at NLO accuracy, both at fixed order and including the matching with parton showers. Relevant differential distributions are studied for two values of the charged boson mass and a thorough comparison is performed between predictions obtained in the four- and five-flavour schemes. We show that the agreement between the two schemes is improved by NLO(+PS) corrections for observables inclusive in the degrees of freedom of bottom quarks. We argue that the four-flavour scheme leads to more reliable predictions, thanks to its accurate description of the bottom-quark kinematics and its small dependence on the Monte Carlos, which in turn is rather large in the five-flavour scheme. A detailed set of recommendations for the simulation of this process in experimental analyses at the LHC is provided.


Introduction
Charged Higgs bosons appear in several extensions of the scalar sector of the Standard Model (SM). In particular, as the SM does not include any elementary charged scalar particle, the observation of a charged Higgs boson would necessarily point to the presence of new physics.
In this paper we focus on one of the simplest extensions of the SM featuring a charged scalar, namely the two-Higgs-doublet model (2HDM). Within this class of models, the existence of five physical Higgs bosons, including two (mass-degenerate) charged particles H ± , is foreseen. Imposing flavour conservation, there are four different ways to couple the SM fermions to the two Higgs doublets. Each of these four ways of assigning the couplings gives rise to a different phenomenology for the charged Higgs boson. Here we consider the so-called type-II 2HDM, in which one doublet couples to up-type quarks and the other to the down-type quarks and the charged leptons. The Minimal Supersymmetric Standard Model (MSSM), up to SUSY corrections, belongs to the type-II 2HDM category.
Light charged Higgs scenarios are classified by Higgs boson masses smaller than the mass of the top quark (typically m H ± 160 GeV), where the top decay to a charged scalar and a bottom quark is allowed. One refers to heavy charged Higgs bosons, on the other hand, for masses larger than the top-quark mass (typically m H ± 180 GeV). In this case, the dominant H ± production channel at the Large Hadron Collider (LHC) is the associated production with a top quark/antiquark and a (possibly low transverse momentum) bottom antiquark/quark. In the intermediate-mass range (160 m H ± 180 GeV) width effects become important and the full amplitude for pp → H ± W ∓ bb (including non-resonant JHEP10(2015)145 contributions) must be taken into account for a proper description. The intermediate-mass range has not been studied so far at the LHC Run I.
Searches at LEP have set a limit m H ± > 80 GeV on the mass of a charged Higgs boson for a type-II 2HDM [1]. The Tevatron limits [2,3] have been superseded by results of the LHC experiments. Recent ATLAS results [4] for a type-II 2HDM based on 19.5 fb −1 of pp collisions at 8 TeV exclude BR(t → bH + )·BR(H + → τ + ν τ ) larger than (0.23−1.3)% in the low-mass region (80 GeV < m H ± < 160 GeV). For the first time limits on tH + production times BR(H + → τ + ν τ ) are provided for 180 GeV < m H ± < 1000 GeV. Rates above 0.76 pb at low masses and 4.5 fb at large masses are excluded. A preliminary note by CMS [5], based on the analysis of 19.7 fb −1 of pp collisions at 8 TeV in the same search channel, sets an upper limit on BR(t → bH + ) · BR(H + → τ + ν τ ) between (0.16 − 1.2)% in the low-mass region (80 GeV < m H ± < 160 GeV). This limit supersedes the one published in ref. [6]. In the high-mass region (180 GeV < m H ± < 600 GeV), limits on tH + production times BR(H + → τ + ν τ ) are set between 0.38 pb and 26 fb. Finally, CMS has also published a preliminary note [7] on a direct search for a heavy charged Higgs which decays in both the H + → tb and the H + → τ + ν τ channels. 1 In this work, we consider heavy charged Higgs boson production at hadron colliders and leave the intermediate-mass range to future studies [9]. In particular, we focus only on the production of a negatively charged scalar since the results are identical for a positively charged scalar. As for any process involving bottom quarks at the matrix-element level, two viable schemes exist to compute the production cross section of a heavy charged Higgs boson. These are usually dubbed as four-and five-flavour schemes. In the four-flavour scheme (4FS) the bottom quark mass is considered as a hard scale of the process. Therefore, bottom quarks do not contribute to the proton wavefunction and can only be generated as massive final states at the level of the short-distance cross section, entailing that b-tagged observables receive contributions starting at leading order (LO). In practice, the theory which is used in such a calculation is an effective theory with four light quarks, where the massive bottom quark is decoupled and enters neither the renormalisation group equation for the running of the strong coupling constant nor the evolution of the parton distribution functions (PDFs). The LO partonic processes in the case at hand are gg → H −b t and qq → H −b t . (1.1) Next-to-leading order (NLO) calculations for the total cross sections in this scheme have been presented in refs. [10,11]. Conversely, in five-flavour scheme (5FS), the bottom quark mass is considered to be much smaller than the hard scales involved in the process. The simplest definition of the 5FS -that suits particularly well perturbative computations -is to strictly set m b = 0 in the short-distance cross section. Consequently, bottom quarks are treated on the same footing as all other massless partons. The only difference is the presence of a threshold in the bottom-quark PDF and the initial condition of the bottom quark evolution being of perturbative nature. The use of b-PDFs comes along with the approximation that, at JHEP10(2015)145 leading order, the massless b quark has a small transverse momentum. In this scheme, the leading logarithms associated to the initial state collinear splitting are resummed to all orders in perturbation theory by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution of the bottom densities. The LO partonic process is given by (1.2) Next-to-leading order predictions for heavy charged Higgs boson production in the 5FS, possibly including the matching to parton-shower Monte Carlos, were studied in refs. [12][13][14][15][16][17][18]. Electroweak corrections [19,20] and soft gluon resummation effects [16,21,22] have also been included in recent works. 2 The leading-order diagrams in the 4FS and 5FS are displayed in figure 1. The comparison between the two schemes at the level of total cross section has been performed by several groups, see e.g. ref. [11] and references therein. In a more recent study [25] a thorough combination of all sources of theoretical uncertainties is performed, state-of-the-art PDF sets are used, the new scale-setting procedure proposed in refs. [26,27] is adopted and a Santander-matched prediction [28] of the four-and five-flavour scheme computations is provided. Furthermore, a comparison of five-and six-flavour schemes, assessing the effect of the inclusion of top PDFs, has been recently performed at the level of total cross sections [29]. In summary, these studies establish a generally good agreement of total cross section predictions in different flavour schemes, once a judicious scale choice is made.
Studies at differential level of heavy charged Higgs production in four-and five-flavour schemes are very limited, if not absent. NLO plus Parton Shower (PS) accurate predictions are available only in the 5FS for both MC@NLO [17] and POWHEG [18], and their comparison to 4FS results at the level of differential distributions has never been performed to date. Such a differential study is certainly more involved than the inclusive one, mainly for two reasons. First of all, there is a very large number of possibly relevant observables at the differential level. Second, differences between the two schemes are generally larger than in the inclusive case. As a consequence, an even more careful assessment of the related uncertainties is necessary for distributions. Indeed, observables inclusive in the degrees of freedom of the bottom quarks, such as total cross sections, turn out to be quite similar in

JHEP10(2015)145
four-and five-flavour schemes. On the one hand power-suppressed terms are small for this class of observables. On the other hand collinear logarithms are phase-space suppressed [26] and therefore moderate, unless a very heavy charged Higgs is produced.
As far as exclusive observables are concerned, it is vital to investigate 4FS and 5FS predictions and assess their differences on a case-by-case basis in an unbiased manner, which is one of the main goals of this paper. The generic rationale is the following: if power-suppressed terms are relevant, the 4FS provides a more reliable prediction, while if collinear logarithms are large, the reorganisation of the perturbative series in the 5FS improves the stability of the perturbative expansion.
The comparison of the two schemes is further complicated by the fact that observables related to light and b-flavoured jets are associated with different perturbative orders: while in charged Higgs production one final-state b quark (not accounting for top decays) is already present at LO in the 4FS, a final-state b quark enters the 5FS computation only at NLO. Furthermore, tagging b quarks in a 5FS fixed-order computation necessarily leads to unphysical results, since they can be considered only as part of jets to retain infra-red cancellations. This issue clearly does not affect the 4FS, in which observables related to the b quarks are regulated by the b mass. Nevertheless, when considering b quarks at very large transverse momenta, b-jet observables have to be preferred anyway, otherwise large logarithms log(p T /m b ) would spoil the perturbative convergence.
Matching to parton showers improves most notably the 5FS predictions, since the parton backward evolution of the initial-state bottom quarks turns them into massive final states and further generates b-flavoured hadrons which renders b-tagging realistic. One should bear in mind, however, that the details and the actual implementation of such backward evolution are highly non-trivial -being based on DGLAP evolution equations, which are only LL accurate -and therefore turn out to be widely Monte Carlo dependent (see section 3.3 of ref. [30], and section 3.3 below). Moreover, the necessity of reshuffling the massless into massive bottom quarks may have significant effects on the kinematics of finalstate B hadrons. In the 4FS, the shower primarily improves the description of Sudakovsuppressed small-p T radiation, by resumming leading collinear logarithms to all orders. In both schemes the PS introduces additional power-suppressed terms in the soft region.
In this paper, we present for the first time NLO+PS accurate 4FS predictions and thoroughly compare them to the 5FS ones. Our primary aim is to acquire a detailed understanding and assess which scheme is better suited for the simulation of the charged Higgs production signal, whose search represents a central part of the physics program in (and beyond) Run II at the LHC. While our results are computed for a type-II 2HDM, they are presented in a model-independent way, which makes them directly applicable also to other 2HDM types and allows for the inclusion of the dominant SUSY contributions, e.g. in the MSSM. As will be shown, our analysis confirms and extends the conclusions that have been drawn in previous studies of differential distributions in the four-and fiveflavour schemes, in the context of Higgs production in association with bottom quarks [31], and Higgs and single-top associated production [32]. Similar analyses were performed for single top production [33,34]. Moreover, given that the mass of the charged Higgs boson is still allowed to take almost arbitrarily large values, this process is perfectly suited for a JHEP10(2015)145 case-study of different flavour schemes in QCD in a broad range of scales.
The paper is structured as follows: in section 2 an outline of the calculation is given and details of the implementation of the 4FS calculation at NLO+PS accuracy are provided. Results are displayed in section 3, in which the features of the 4FS calculation are described, including the size of the interference term that appears in this scheme, and the 4FS and 5FS distributions are compared. We conclude in section 4 with our final recommendations for experimental analyses.
2 Outline of the calculation Our computation takes advantage of a full chain of automatic tools developed to study the phenomenology of new physics models at NLO QCD accuracy in the Mad-Graph5 aMC@NLO [35] framework.

Framework
Besides the usual tree-level Feynman rules, some extra ingredients have to be provided to the matrix-element generator in order to obtain the code for the simulation of a new physics process at NLO. These extra ingredients are the Ultra-Violet (UV) renormalisation counterterms and a sub-class of the rational terms that enter the reduction of the virtual matrix elements, the so-called R 2 terms [36]. UV and R 2 counter-terms can be computed starting from the model Lagrangian for any renormalisable theory via the NLOCT package [37], based on FeynRules [38] and FeynArts [39]. The tree-level and NLO Feynman rules are exported in the Universal Feynrules Output (UFO) [40] format, as a Python module which can be loaded by matrix-element generators, such as MadGraph5 aMC@NLO. Finally, the UFO information is translated into helicity routines [41] by ALOHA [42].
We remind the reader that MadGraph5 aMC@NLO is a meta-code that automatically generates the code for simulating any process at NLO(+PS) accuracy. Mad-Graph5 aMC@NLO adopts the FKS method [43,44] for the subtraction of the singularities of the real-emission matrix elements, as automated in the MadFKS module [45]; virtual matrix elements are computed via the MadLoop module [46], which relies on the OPP [47] and Tensor Integral Reduction (TIR) [48,49] methods, as implemented in CutTools [50] and IREGI [51], supplemented by an in-house implementation of the OpenLoops procedure [52]; finally, the generation of hard events and their matching to parton-shower simulations is performedà la MC@NLO [53]. The matching to Her-wig6 [54], Pythia6 [55] (ordered in virtuality or in transverse momentum, the latter only for processes with no light partons in the final state), Herwig++ [56] and Pythia8 [57] is available. As a consequence, the only inputs needed are the implementation of the model in FeynRules and the inclusion of the running of the bottom quark mass, which follows from our renormalisation-scheme choice. More details are given in the next section. Both the 4FS and the 5FS computations have been performed in the framework specified above, with the additional advantage of ensuring the consistency of all settings and input parameters.
To guarantee the validity of our analysis, we have performed a detailed comparison of the inclusive cross sections obtained with MadGraph5 aMC@NLO against previous JHEP10(2015)145 results published in literature. The calculation in the 4FS has been compared to the private code of ref. [11] at the level of total rates, and a good agreement within the numerical uncertainty of the reference code has been found. Furthermore, with the bottom quark Yukawa renormalised on-shell and with its mass set the value of the top pole mass, we reproduce the ttH total cross section in the SM at the few per-mille level.

Implementation
We have used the implementation of the generic 2HDM in FeynRules detailed in ref. [37]. This model has been converted into a type-II 2HDM by adding β as an external parameter and by restricting the Yukawa couplings accordingly. If top and bottom quarks are assumed to be the only massive fermions, the only non-zero entries of the Yukawa coupling matrices to the doublet without vacuum expectation value in the Higgs basis for the type-II 2HDM are given by where m y t/b are the Yukawa masses of the top and bottom quark. The parameter tan β = v 2 /v 1 is the ratio of the vacuum expectation values v 1 and v 2 of the two Higgs doublets, With those restrictions, the H − tb vertex is given by where P R/L = (1 ± γ 5 )/2 are the chirality projectors and y t/b ≡ √ 2 v are the corresponding SM Yukawa couplings. We strictly separate the Yukawa masses that are used in the computation of the couplings between the fermions and the scalars from the kinematic masses that are used everywhere else and set to the on-shell mass. 3 This distinction allows us to keep a non-vanishing bottom Yukawa in the five-flavour scheme as the leading term in the small m b expansion [14,15]. Furthermore, it allows us to choose different renormalisation schemes for the bottom quark mass in the matrix element and in the Yukawa coupling.
The model R 2 and UV vertices required for NLO computations in Mad-Graph5 aMC@NLO have been computed using NLOCT [37]. The masses and the wave functions are renormalised in the on-shell scheme to avoid the computation of loops on external legs. The strong coupling constant is renormalised in the MS scheme with the contribution of massive quarks subtracted from the gluon self-energy at zero-momentum transfer. Therefore, only the massless modes affect the running of α s . The renormalisation of the masses in principle fixes the renormalisation of the top and bottom Yukawa since

JHEP10(2015)145
in the on-shell scheme. This is the default renormalisation used in NLOCT and it would ensure that nothing but the strong coupling constant depends on the renormalisation scale. The top mass and Yukawa are always renormalised in this way throughout this paper. Therefore its Yukawa mass is set equal to the pole mass. On the contrary, the bottom quark Yukawa has been renormalised in the MS scheme, i. e.
This scheme choice has the advantage of resumming potentially large logarithms log(µ R /m b ) (with µ R ∼ m H ± ) to all orders. The bottom Yukawa mass is set to the value of the running MS mass at the renormalisation scale. Besides the modifications at the level of the UFO model, also the code written by MadGraph5 aMC@NLO had to be changed in order to account for the additional scale dependence introduced by the b-quark Yukawa, in particular for what concerns the on-the-fly evaluation of scale uncertainties obtained via reweighting [59]. This has been done in an analogous way as for bottom-associated Higgs production [31], by splitting the cross section in parts that factorise different powers of y b , i. e. y 2 b , y b y t and y 2 t .

Results
In this section, we present four-flavour scheme predictions of charged Higgs boson production at NLO matched to parton showers. This calculation has never been performed before in the literature. Several differential distributions that are reconstructed from the final state particles in tH −b production are studied. We investigate the role of the shower scale in this process, discuss the impact of the y b y t interference term and compare our reference predictions at (N)LO+PS to the f(N)LO results. For matched predictions, both Herwig++ and Pythia8 are employed. We conclude this chapter with a comprehensive comparison of 4FS and 5FS distributions, in which the effects of higher order corrections, the impact of the choice of the shower scale and the dependence of each scheme on the different Monte Carlos are analysed.

Settings
We present results for charged Higgs boson production at the LHC Run II ( √ S had = 13 TeV) by considering two scenarios: a lighter (m H ± = 200 GeV) and a heavier (m H ± = 600 GeV) charged Higgs boson. For simplicity, we set tan β = 8 throughout this paper. At this value, y 2 b and y 2 t terms are of similar size and the relative contribution of the y b y t term to the total cross section is close to its maximum. Results for any other value of tan β can be obtained by a trivial overall rescaling of the individual contributions according to their Yukawa couplings (y b by tan β, y t by 1/ tan β). Therefore, we preserve the generality of our results by studying the y 2 b , y 2 t and y b y t contributions separately. Our results are further applicable to supersymmetric (SUSY) models, such as the Minimal Supersymmetric SM (MSSM), in which case additional higher-order contributions associated with the virtual exchange of supersymmetric particles have to be included. Of JHEP10(2015)145 particular relevance are corrections that modify the tree-level relation between the bottom quark mass and its Yukawa coupling. These corrections are enhanced at large tan β. The leading contribution to SUSY corrections in charged Higgs boson production is taken into account via a modification of the Yukawa coupling that resums tan β-enhanced terms to all orders, see refs. [11,[60][61][62][63][64][65][66][67][68][69][70]. For some MSSM searches, they may significantly change the interpretation of the discovery/exclusion limits [65,66]. Splitting our predictions in terms that factorise the different powers of the Yukawas, therefore, allows us to easily implement this class of SUSY corrections via a rescaling of the Yukawa couplings. The remaining SUSY-QCD effects are marginal at large tan β, but can reach up to O(10%) at small tan β [11]. Any supersymmetric correction is, however, beyond the scope of this paper and will not be discussed any further.
We show results obtained with the NNPDF2.3 set [71] at NLO and the NNPDF3.0 [72] set at LO. To obtain consistent predictions, parton distribution functions (PDFs) computed in the proper flavour number scheme are used: we interface our NLO (LO) calculation with the NNPDF2.3 (NNPDF3.0) with n f = 4 and n f = 5 active flavours for the 4FS and 5FS respectively. The mismatch between the PDF sets used in the LO and NLO computations is due to the absence of a public set of non-QED LO PDFs in the NNPDF2.3 family. This does not affect the accuracy of our results, given that the LO PDF sets exhibit a theoretical uncertainty which is larger than the difference between the two NNPDF families. The strong coupling constant is consistent with α s (M Z ) = 0.118 for the 5FS NLO parton densities and α s (M Z ) = 0.1226 for the 4FS NLO ones. 4 The heavy quark pole masses are set to Finally, our central renormalisation and factorisation scales µ R , µ F are set to where the index i runs over all final state particles (the top quark, the charged Higgs boson and possibly the extra b quark and/or light parton) of the hard process. For vanishing transverse momenta of the external particles, our scale choice corresponds to the factorisation scale set in the 4FS calculation of refs. [11,25]. In the following, scale uncertainties are obtained by varying µ F and µ R independently by a factor of two around their central values, given in eq. (3.3). We have checked that, particularly for our reference 4FS NLO+PS prediction, the dependence of the distributions on the shower scale µ sh , when varied by a

JHEP10(2015)145
factor in the range [1/ √ 2, √ 2], is rather mild and significantly smaller than uncertainties associated with the renormalisation and factorisation scales; we therefore refrain from including uncertainties associated with µ sh in what follows. Furthermore, we will not discuss any PDF systematics. 5 Jets are reconstructed via the anti-k T algorithm [73], as implemented in FastJet [74,75], with a distance parameter ∆R = 0.4 and subject to the conditions For fixed-order computations jets are clustered from partonic final states, while in simulations matched to parton showers jets are made up of hadrons; b jets are defined to contain at least one b quark (at fixed order) or B hadron (in matched simulations). In our simulations we keep the charged Higgs boson stable, while we decay the top quark leptonically (although the leptons from the decay will not affect any observable we consider) in order to keep as much control as possible on the origin of QCD radiation. The task to decay the top quark is performed by the parton shower for (N)LO+PS runs, while at fixed order we simulate the decay t → bW in an isotropic way (in the t rest frame) at the analysis level. 6 No simulation of the underlying event is performed by the parton shower.
Let us conclude this section by addressing one further point, which is crucial when processes with final-state b quarks are matched to parton showers: the choice of the shower starting scale µ sh . Such processes are known to prefer much lower values of the renormalisation and factorisation scales than the one naively identified as the hard scale of the process (ŝ). In fact, the shower starting scale and the factorisation scale emerge both from the same concept, namely the separation of soft and hard physics. Furthermore, it has been argued in ref. [31] for the associated production of a neutral Higgs boson with bottom quarks that the shower starting scale (limiting the hardest emission that the shower can generate) should be set at similar values, i. e. well belowŝ. Following the arguments made in ref. [31], we check their validity in the case of charged Higgs boson production. We shall stress at this point that the following discussion applies both to our reference scenarios with m H − = 200 GeV and m H − = 600 GeV, although we refrain from showing explicit results for the latter.
MadGraph5 aMC@NLO assigns a dynamical shower scale chosen from a distribution in the range 7 where F is a parameter that drives the bounds of the distribution, and whose default value is F = 1. With such a default setting the effective value of µ sh , namely the maximum of the µ sh distribution (which for simplicity we will refer to as just µ sh in the following), is indeed 5 Note that scale variations due to µF and µR as well as PDF uncertainties are computed at no extra CPU cost using the reweighting procedure of ref. [59]. 6 Such an approach neglects spin-correlation in the decay of the top quark. However, within the Mad-  term, orange dotted-solid for the y 2 t term), and at NLO+PS with F = 1 (green dots for the y 2 b term, orange dots for the y 2 t term) and F = 4 (green solid for the y 2 b term, orange solid for the y 2 t term). We show predictions matched with Pythia8 (left) and Herwig++ (right). The insets show the ratio of the curves in the main frame over the fNLO prediction, for both the y 2 b and the y 2 t terms.
much larger than µ F,R . Furthermore, considering the transverse momentum distribution of the Born-level "system" (p T (sys)), 8 which is maximally sensitive to the interplay between the fixed-order prediction and the shower, the NLO+PS distribution (in particular in the 4FS) does not match the fixed-order NLO (fNLO) one at large p T for F = 1. This can be deduced from figure 2, when comparing the crosses (NLO+PS for F = 1) to the solid curves overlayed with points (fNLO). On the contrary, we observe a clearly improved highp T matching of the NLO+PS results to the fixed-order ones by choosing a reduced shower scale corresponding to F = 4 (solid curves). 9 Indeed, such a choice brings µ sh much closer to the value of the renormalisation and factorisation scales. We have also checked that the agreement among Pythia8 and Herwig++ improves (although often only marginally) when differential observables in the 4FS are computed with F = 4.
In conclusion, although for this process we do not reproduce all results of ref. [31] with the same significance, we still find sufficient evidence that F = 4 is favourable in many 8 Note that the Born-level system is unambiguously defined only in a fixed-order calculation, being in our case the charged Higgs accompanied by the final state top and bottom quark. At NLO+PS we define it to include the hardest B hadron (instead of the bottom quark), which does not originate from the top decay; in this case, MC-truth is used. 9 Our focus here is on the 4FS prediction. However, similar conclusions, if less stringent, can be drawn from the corresponding plots in the 5FS. respects and make it our default choice. In section 3.3, we shall further study the impact of this choice when comparing the 4FS and 5FS results: by setting F = 4 an improved agreement between the two schemes at the level of shapes is observed.

Four-flavour scheme results
We now turn to our phenomenological results for charged Higgs boson production. Let us first consider state-of-the-art 4FS predictions, which, as will be shown, constitute the most reliable differential results for observables exclusive in the degrees of freedom of finalstate bottom quarks. We split this section into two parts: in section 3.2.1 we limit our study to the dominant y 2 b and y 2 t contributions, while the y b y t contribution is considered in section 3.2.2.

y 2
b and y 2 t contributions at NLO+PS We begin our analysis by studying total rates for the production of charged Higgs bosons with a mass of 200 GeV and 600 GeV in table 1. We consider three possibilities: the fully inclusive case, the case in which we require at least one b jet, and the one in which two or more b jets are tagged. All results are given at both LO and NLO accuracy. The cross sections in which one or two b jets are required depend on the approximation and Monte

JHEP10(2015)145
Carlo under consideration. We thus report separately results obtained at fixed order, with Pythia8 and with Herwig++. Any quoted uncertainty is due to scale variation, evaluated as detailed in section 3.1; they are indicated only at fixed order and for results matched with Pythia8, since they show little dependence on the specific Monte Carlo. Results for y 2 b and y 2 t terms are presented separately. Let us summarize the conclusions to be drawn from table 1 as follows: • The scale uncertainty of NLO predictions is substantially smaller than that of the LO ones; at NLO the scale uncertainty is larger for the y 2 b than for the y 2 t contribution (∼ 15-20% and ∼ 10-15%, respectively), due to the different renormalisation schemes used for the bottom and top Yukawa couplings.
• Because of our default choice of tan β = 8, y 2 b and y 2 t predictions are of similar size at NLO (only ∼ 15% different); the difference is larger at LO (∼ 30%). As a consequence, the K-factors are generally different between the y 2 b and y 2 t terms; for m H − = 200 GeV, the inclusive y 2 b K-factor is close to 1.2, while for the y 2 t term the NLO corrections have a larger impact, with K ≈ 1.5; for m H − = 600 GeV, NLO corrections are larger in both cases: K ≈ 1.3 and K ≈ 1.6 for the y 2 b and y 2 t terms respectively. The smaller K-factor for the y 2 b term is a consequence of the fact that, by renormalising the bottom Yukawa in the MS scheme, the LO predictions already include a class of higher order corrections.
• Pythia8 and Herwig++ predictions for cross sections with b-jet requirements are consistent with each other and well within quoted uncertainties; this holds true both at the LO and at the NLO. Even the agreement with the fixed-order results is rather good overall and largely within uncertainties. Higher b-jet multiplicities tend to deteriorate the agreement, although only to a small extent.
• The effect of the cuts on the number of b jets is quite moderate: given that there is generally a hard b jet from the top quark decay, the inclusive rate is reduced by only ∼ 10%, when at least one b jet is required. On the other hand, requiring a second b jet reduces the cross section by a factor of 4-5. Scale uncertainties slightly decrease at NLO for cross sections within cuts.
We now turn to differential observables in the 4FS. We only discuss results at fixed order and matched with Pythia8. Any differences to the matched Herwig++ predictions will be explicitly mentioned below. Besides, the Monte Carlo dependence will be investigated in more detail when comparing distributions in the two flavour schemes in section 3.3.
The figures throughout this section are organised according to the following pattern: the main frame displays the relevant predictions in absolute size as cross section per bin 10 for the y 2 b and y 2 t terms. For the sake of readability, they have been separated by reducing the y 2 t curve by a factor of ten. The fixed-order results are shown using boxes, while JHEP10(2015)145  normalised to the corresponding NLO+PS result. The difference between the two insets are the uncertainty bands, which in the first one are displayed for the LO predictions, while in the second for the NLO ones. The third and fourth insets are identical to first two, but for the y 2 t contribution. Let us start in figure 3 with the transverse momentum spectrum of a 200 GeV Higgs (top left panel) and the associated top quark (top right panel). The two plots are in fact quite similar and can be discussed simultaneously: the agreement between the fixed-order and the PS-matched simulations is close-to-excellent in all cases. Matching to the PS has only minor effects on the shape of the distributions, notable only at small transverse momenta. In other words, the resummation effects of the shower are extremely small, as can be expected from observables that are practically not affected by any large Sudakov logarithm. The NLO corrections mostly affect the normalisation and reduce theoretical uncertainties, reflecting the numbers quoted in table 1. We stress here two general features regarding the relation between the y 2 b and y 2 t terms which shall hold true for all subsequent observables under consideration: first, the MS renormalisation makes the uncertainty associated with the y 2 b curves larger, due to the variation of µ R in y b ∼m b (µ R ); second, the ratio of the y 2 b and y 2 t contributions at NLO is generally flat. More details will be given in section 3.3.
For a heavier charged Higgs boson with mass m H − = 600 GeV, the lower panels of figure 3 display a similar behaviour as in the m H − = 200 GeV case. Larger effects due to the PS are visible in the Higgs p T spectrum at low transverse momenta.
We continue our presentation of the 4FS results with the transverse momentum spectra JHEP10(2015)145 of the hardest (p T (b 1 )) and second-hardest (p T (b 2 )) b jet in figure 4; see section 3.1 for our jet (and b jet) definition. In this case, we limit our discussion to m H − = 200 GeV results, since, apart from a naturally smaller cross section, the relative behaviour of the different curves is extremely similar for a heavier charged scalar. Although the p T (b 1 ) and p T (b 2 ) distributions in the main frame develop a quite different behaviour in terms of their shapes, the ratios in the insets exhibit comparable features: in all cases, the tail of the spectra is driven by the order relevant to the simulation and NLO corrections slightly soften the spectra. In other words, the fixed-order results agree rather well with the corresponding Pythia8 ones in the tail. On the other hand, close to threshold (p T = 25 GeV), where resummation effects are enhanced, non-showered and showered results exhibit sizeable differences, in particular for the hardest b jet.
Turning now to somewhat related observables in figure 5-the transverse momentum distributions of the hardest and second-hardest B hadron-, one may expect rather similar features to the b-jet transverse momentum spectra. On the contrary, their pattern is actually very different; the salient feature being that showered results are vastly softer than the fixed-order ones and a substantial shape distortion due to the matching with parton showers is observed. In fact, even the peak of the p T (B 1 ) distribution is moved by ∼ 25 GeV towards the left by the shower. However, one should bear in mind that we compare bottom quarks at parton level for the f(N)LO predictions with B hadrons at (N)LO+PS. The observed differences unavoidably lead to the conclusion that fragmentation effects become significant for such exclusive observables. Otherwise, the pattern of the Pythia8 results is JHEP10(2015)145 very much reminiscent of b-jet spectra, displaying a slightly harder LO+PS shape than at NLO+PS. Generally speaking, the relative behaviour of the y 2 b and y 2 t curves is pretty much alike, including the peculiar increase of the f(N)LO cross section towards vanishing p T (B 2 ). Again, we refrain from showing explicit results for a m H − = 600 GeV charged Higgs boson, since the pattern of the various curves turns out to be very similar to the m H − = 200 GeV case; the only difference to be pointed out is a slightly reduced gap between showered and fixed-order results for m H − = 600 GeV.
We investigated a vast number of differential observables, the majority of which follows the same pattern as illustrated in figures 3 and 4: the NLO corrections are rather flat and lie within the LO uncertainty bands, shower effects are moderate and become more substantial the more exclusive the observable is with respect to the bottom-quark degrees of freedom. Based on our findings, we do not expect larger effects for b-jet observables than for distributions relevant to B hadrons. Therefore it is worth to consider the invariant mass M (B 1 , B 2 ) of and the distance ∆R(B 1 , B 2 ) in the η − φ plane between the two hardest B hadrons, which are displayed in the left and right panels of figure 6 respectively. The reader should keep in mind that in the fixed-order cases the two hardest b quarks are used instead. Since there are no salient differences for larger Higgs masses we only discuss the m H − = 200 GeV results.
Although the invariant mass is quite a different observable than the transverse momentum of the hardest B hadron, our findings are actually rather similar, but less pronounced: the shower substantially affects the distributions, by causing a significant softening whose size exceeds the uncertainty bands of the fixed-order calculations. This effect reflects the JHEP10(2015)145 loss of energy of the b quarks due to their fragmentation into B hadrons. These observations are to a good extend independent of the considered contribution (y 2 b or y 2 t ). In contrast, the effect of the shower in the ∆R(B 1 , B 2 ) distribution is smaller and becomes relevant only at large separations (∆R ≥ 4.5). We point out that such large separations are at or beyond the coverage edge of the b tracking system of LHC experiments, and that if two b jets are explicitly required these differences are reduced. Effects due to the parton shower are also visible at small separations, where secondary g → bb splittings can be important. In this case, the two B hadrons are likely to be clustered in the same b jet. The salient differences between y 2 b and y 2 t contributions are related to the behaviour of the respective LO curves, which again depends on their different relative normalisations; their shapes, on the contrary, are quite similar.
Let us now conclude our 4FS study by considering jet rates, displayed in figure 7. In the left panel, we show jet multiplicities without requiring any b tagging, while in the right panel we show b-jet multiplicities. First of all, we recall that in our setup the top quark decays (leptonically) both in the shower and at fixed order. For this reason, up to two/three jets and up to two b jets can appear at fLO/fNLO. Looking at the m H − = 200 GeV results first (upper plots), we can appreciate the effects of the NLO corrections and of the matching with parton showers. NLO effects have the largest impact in the two-jet bin, where the cross section is increased by 35% (65%) for the y 2 b (y 2 t ) term. Their effect is minor for lower multiplicities: almost no correction occurs for the y 2 b term, while the y 2 t one is increased by 20%. If in turn we consider the effects of the parton shower and compare the fNLO histograms to the NLO+PS ones, we infer that the main effect is to populate higher multiplicities. Therefore, the zero-and one-jet rates are reduced as a consequence of unitarity. Such a reduction is quite important, as the NLO+PY8 prediction falls outside the fNLO uncertainty band, particularly in the one-jet bin. The two-jet bin is left almost unchanged by the shower. We also find that shower effects have a much larger impact at LO, which reflects the large uncertainties associated with the LO computations, particularly at fixed order.
The distribution of events with respect to the number of b jets is displayed in the upper right panel of figure 7. In this case NLO corrections have the largest effect on the zero-and one-b jet bins. Since the majority of events lies in the one-b jet bin, the shower moves events from this bin into the higher and lower multiplicities, although the effect is in general moderate and within the fNLO uncertainties. Only in the zero-b jet bin, whose rate is however suppressed, the differences between fNLO and NLO+PS results are larger than in the flavour-unspecific case and the uncertainties barely overlap. Overall, the NLO predictions are reasonably close to each other, since higher multiplicities (beyond two b jets) are phase-space suppressed in the NLO+PS simulations.
For both the jet and b-jet rates, the fraction of events with jet multiplicities beyond the ones already present at the hard-matrix element level (more than three jets and more than two b jets) is in good agreement between the LO+PS and NLO+PS predictions, the LO ones being slightly enhanced though. At this point, we remark that these multiplicities are utterly Monte Carlo dependent and a substantial disagreement between Pythia8 and Herwig++ is found for these bins. For the lower multiplicities, however, their agreement JHEP10(2015)145  normalisation, the distribution of events at LO appears shifted towards low multiplicities as compared to LO+PS, without any overlap of the corresponding uncertainties in the first two bins. However, given our findings for the lower Higgs mass case, this is expected: the shower shifts events from lower towards higher jet multiplicities; this is enhanced for m H − = 600 GeV due to a generally increased hardness of the process. Indeed, the two-jet bin has the largest rate, and the three-and four-jet bins are less suppressed than in the lighter Higgs case. NLO corrections slightly improve the agreement of showered and fixedorder results, albeit fNLO and NLO+PS still fall outside the respective uncertainties in the zero-and one-jet bins.
In the case of b jets, on the other hand, the features of the relative curves reflect those discussed for the lighter Higgs and no further comments are needed.

The y b y t contribution
As we mentioned before, in the 5FS the NLO cross section receives contributions either proportional to y 2 b or to y 2 t . No y b y t term appears, given that it would come from the interference of left-handed with right-handed massless bottom quarks. If in turn b quarks are massive, as in the 4FS, the y b y t term does not vanish any longer, and it is proportional to m 2 b /Q 2 , where Q is some hard scale of the process. So far, we have limited our 4FS analysis to the y 2 b and y 2 t contributions, assuming the y b y t one to be suppressed. In this section, we show that this is indeed the case.
To this purpose, we consider the total cross section for charged Higgs production in the 4FS at LO and NLO, and plot the relative contribution −σ y b yt /σ all as a function of tan β in figure 8, with σ all being the sum of all terms. The results are shown for m H − = 200, 600 GeV. The minus sign takes into account the fact that the y b y t term is negative. We stress that the y b y t contribution is independent of tan β. As can be inferred from the plots, the relative size of the y b y t term is below 5% for m H − = 200 GeV, and 0.5% for m H − = 600 GeV. The relative contribution to the cross section proportional to y b y t is JHEP10(2015)145 Let us further investigate the potential impact of the inclusion of the y b y t term on some differential observables, for such a value of tan β. In particular, we look at the transverse momentum of the Higgs, the top and the two hardest B hadrons for m H − = 200 GeV, displayed in figure 9. From these plots we notice that the effect of the y b y t term is peaked at low scales, by reaching at most 6 − 7% of the full cross section, and is almost the same 11 The difference between the LO and NLO values is due to the different perturbative order in the running of y b . at LO and NLO. We stress again that these numbers have been computed for the value of tan β for which the relative y b y t contribution is maximal: for larger (smaller) values of tan β, this contribution is suppressed by a factor 1/ tan 2 β (tan 2 β) with respect to the y 2 t (y 2 b ) contribution and further reduced for heavier charged Higgs bosons. The typical scale uncertainties at NLO (∼ 10 − 15%) justify our choice to neglect the y b y t contribution in the current analysis. A viable alternative would be to include the relative contribution of the y b y t term only at LO, which was shown to be very similar to the NLO one.

Four-and five-flavour scheme comparison
We turn now to investigate how predictions obtained in the four-and five-flavours schemes compare. The two schemes are actually identical up to b-mass power suppressed terms when computed to all orders in perturbation theory, but the way of ordering the perturbative series is different. As a consequence, the results in the two schemes may be different at any finite order, while the inclusion of higher orders necessarily brings the predictions in the two schemes closer to each other. We start by quantifying how the inclusion of NLO corrections improves their mutual agreement. In figures 10-12 we show, for some relevant observables, the LO and NLO predictions (matched with Pythia8) in the two schemes. All figures have the same pattern: a main frame with the absolute predictions in the 5FS JHEP10(2015)145  frame over the 5FS NLO prediction, for the y 2 b and y 2 t contributions respectively. For the NLO predictions, a band indicating the scale uncertainty 12 is attached to the curves. In the third inset, the four differential NLO/LO K-factors (y 2 b and y 2 t for 4FS and 5FS) are displayed.
A general observation is that, as expected, the NLO predictions in the two schemes are much closer to each other than the LO ones, in particular as far as shapes are concerned. Differences in the overall normalisation reflect the differences in the total cross section, which were already discussed in ref. [25], while in this comparison we are mostly interested in the shapes. In figure 10 we observe that for the transverse momentum of the top quark and the Higgs boson the difference between the two schemes can be compensated by a simple overall rescaling of the total rates (σ 4FS tot /σ 5FS tot 0.7) at NLO, while LO predictions in the two schemes have considerably different shapes. The same level of agreement should be found also for observables related to the (leptonic) decay products of the top quark and the Higgs. Let us recall that in our simulation we do not decay the Higgs boson, but we decay leptonically the top quark. The b quark from the top decay mostly ends up in the hardest b jet. This explains why the p T spectrum of the hardest b jet (left plot in figure 11) displays a flat ratio between the 4FS and 5FS at NLO, up to ∼ 120 GeV. Above that value, secondary g → bb splittings from hard gluons become more relevant, which is also reflected in the growth of the 5FS uncertainty band and K-factor. A similar behaviour has been observed in the case of tH production in the SM [32]. The pseudo-rapidity of the hardest JHEP10(2015)145 b jet (right plot in figure 11) is mostly dominated by the low-p T region, and it therefore also displays a good agreement between 4FS and 5FS shapes at NLO.
Larger differences between the two schemes appear for the second-hardest b jet, which is expected to be poorly described in the 5FS. In particular, its kinematics in the 5FS at LO is determined by the shower, while at NLO it is driven by a tree-level matrix element (therefore being formally only LO accurate). Predictions for the transverse momentum of the second b jet and its pseudo-rapidity are shown in the left and right panels of figure 12. The 5FS develops large K-factors and larger uncertainties, since its LO prediction stems only from the shower evolution. Therefore, the 4FS description has to be preferred for these observables, both because of its better perturbative behaviour and the proper modelling of the final-state b jets.
The effects of the different treatment of the bottom quark in the two schemes is even more visible for the differential observables related to the hardest B hadron (see figure 13). At medium and large p T (B 1 ) and at central η(B 1 ) similar effects as for the hardest b jet are observed. At variance, the 4FS prediction is suppressed with respect to the 5FS one at low p T (B 1 ) and at large η(B 1 ). This is most likely due to mass effects: these kinematical regions correspond to one b quark being collinear to the beam. In the 5FS these configurations are enhanced because of the collinear singularities, while in the 4FS such a singularities are screened by the b-quark mass. Therefore, even after the PS, the 5FS is reminiscent of the collinear enhancement. In the case of the second-hardest B hadron (not shown) these effects are further enhanced.
Let us make a final remark on the inclusion of the NLO corrections. The NLO/LO K-factor is quite different in the two schemes: in the 4FS the K−factor appears much more pronounced for the y 2 t than for the y 2 b term, while in the 5FS it is similar for both contributions. Despite that, a remarkable compensation in shape between the LO differential cross sections and the NLO corrections takes place, such that the 4FS/5FS ratio at NLO is quite similar for the y 2 b and y 2 t terms. All the plots discussed so far are relevant to the lighter Higgs under consideration (m H − = 200 GeV). For a heavier charged Higgs boson (m H − = 600 GeV), the picture does not change significantly. The only thing that may be worth mentioning is the fact that, for the second b jet, the K-factor in the 5FS lies much closer to unity than for the lighter Higgs. Such a behaviour may be due to the increased weight of the initial-state collinear logarithms resummed by the bottom PDFs in the 5FS, which are enhanced at larger masses of the produced particle. Besides, as already pointed out in ref. [26], collinear logarithms become increasingly relevant the larger the fraction of the momentum carried by the initial partons is.
We continue our analysis by investigating how the choice of the shower scale affects the results in the two schemes. We stress once more that our default choice (F = 4) is physically well-motivated by the arguments given in section 3.1. Below we show that this choice also improves the mutual agreement between the NLO predictions in the two schemes. A number of differential distributions in the 4FS and 5FS for F = 1 and  4FS (in blue and orange respectively). Solid curves are used for our reference predictions with F = 4, while dashed curves refer to the default choice in MadGraph5 aMC@NLO (F = 1). The first inset displays the ratio of the curves in the main frame over the corresponding ones (y 2 b or y 2 t ) in the 5FS for F = 4. The second inset shows the ratio of the curves with F = 1 in the 4FS over the corresponding ones in the 5FS. In these two insets we can analyse whether F = 1 or F = 4 yields a flatter 4FS/5FS ratio. In the last inset, the ratio of the normalised y 2 b and y 2 t distributions is plotted, for the 4FS and 5FS and for F = 1, 4. The purpose of this inset is to study whether the y 2 t and y 2 b contributions develop similar shapes.
Overall, we observe a smaller dependence on F in the 5FS than in the 4FS distributions. A similar behaviour was found in the context of Higgs production in association with bottom quarks [31]. In any case, the choice of F is almost irrelevant when considering observables that are not sensitive to the b-quark degrees of freedom. In figure 14, where m H − = 200 GeV, we notice that charged Higgs and top p T distributions are slightly harder in the 4FS (10-15% at large p T ) for F = 1. Moreover, the F = 4 choice yields a remarkably flat 4FS/5FS ratio, much flatter than for F = 1. Consequently, F = 4 improves the agreement between the two schemes.
This improvement becomes even more visible for observables sensitive to the b kine-JHEP10(2015)145   4FS than F = 4. A similar pattern is visible in the case of B hadrons. In this case, however, the general agreement between 4FS and 5FS significantly deteriorates, as it was already observed.
In general, similar conclusions can be drawn for m H − = 600 GeV. The only two distributions that exhibit a different behaviour as compared to the m H − = 200 GeV case are shown in figure 16. The transverse momentum distribution of the heavy charged Higgs boson (left panel) in the 4FS is much more affected by the choice of F , with effects that reach up to 40% in the tail of the distribution. In this case the 4FS and 5FS mutual agreement is significantly improved by the choice F = 4. The transverse momentum distribution of the second-hardest b jet (right panel) is less sensitive to the choice of F than in the lighter-Higgs scenario.
Jet multiplicities are a class of observables that are particularly sensitive to the excess of radiation generated by using F = 1. From figure 17 we conclude that the 4FS generally prefers higher jet multiplicities than the 5FS. A similar behaviour was observed also in ref. [32] and it was explained by considering the different colour structure of the initial state in the two schemes, besides the fact that the process is generally harder in the 4FS. This tendency is slightly reduced -and as a consequence the agreement slightly improves -by the choice of a smaller shower scale (F = 4). We remark that the dependence on the shower scale is even less apparent in the heavy charged Higgs case, and for b-jet multiplicities (not shown).
Finally, the ratio between the y 2 b and y 2 t terms (last insets of figures [14][15][16][17] illustrates that the two contributions give remarkably similar shapes for all observables under con-JHEP10(2015)145 Figure 18. Transverse momentum (left) and pseudo-rapidity (right) distribution of the secondhardest B hadron. NLO curves matched with Pythia8 (solid) and Herwig++ (dashed) are shown, in the 5FS (blue for the y 2 b term, dark-red for the y 2 t one) and 4FS (green for the y 2 b term, light-orange for the y 2 t one). The two insets show the ratios of the histogram in the main frame over the 5FS prediction matched with Pythia8, for the y 2 b and y 2 t terms separately. Scale uncertainty bands are shown around the Pythia8 predictions. sideration. Any difference emerges from the dynamical scale choice at which the bottom Yukawa is computed (m b (µ R ) with µ R = H T /3, see section 3.1). However, such difference is never larger than 10%, with the y 2 b distributions being slightly softer than the y 2 t ones. Finally, we analyse the sensitivity of various observables to the parton shower in the four-and five-flavour schemes. To this purpose, we compare results at NLO matched with the Pythia8 and Herwig++ Monte Carlos. Having verified that the relative behaviour of the two Monte Carlos hardly depends on the specific choice of the charged Higgs mass under consideration, we limit the discussion to the m H − = 200 GeV results.
As already pointed out in the introduction, the parton-shower matching for bottomquark initial states in the 5FS involves some approximations: the initial-state backward evolution is based on leading-log accurate gluon splittings and requires the reshuffling of massless into massive bottom quarks. For these reasons, the 5FS predictions are extremely sensitive to the specific treatment of bottom quarks in a given Monte Carlo. This is most remarkable for b-jet/B-hadron related observables, such as the transverse momentum distribution of the second-hardest B hadron (displayed in the left panel of figure 18). In the 5FS, the Herwig++ prediction displays a significant shape distortion, both at small and large values of the transverse momentum, falling outside the Pythia8 uncertainty bands at small p T . A similar discrepancy is evident in the forward region of the rapidity spectrum, shown in the right panel of figure 18. In this plot, the 5FS prediction is larger for Herwig++ than for Pythia8. For both observables, the 4FS results display a significantly JHEP10(2015)145 smaller Monte Carlo dependence.
An even more spectacular example is provided by the η − φ distance between the two hardest B hadrons. The results are shown in the left panel of figure 19. The Herwig++ prediction in the 5FS features a much higher tail at large separations. This can be traced back to the fact that Herwig++ tends to produce B hadrons much closer to the beam line when simulations are performed in the 5FS. The same behaviour was observed in ref. [31], in which it was pointed out that such effects are however not relevant when realistic cuts on the B hadrons are imposed, for example, when b jets are required. As a matter of fact, we observe a neatly improved agreement among all curves in figure 19 when requiring at least two b jets (right panel).
Looking further at jet (left panel) and b-jet (right panel) multiplicities in figure 20, we observe instead that the two flavour-schemes display a similar Monte Carlo dependence. In the flavour-unspecific case this dependence is quite small for all jet-multiplicity bins. On the other hand, for b jets Herwig++ and Pythia8 are in good agreement up to one b jet in the 5FS and two b jets in the 4FS, that is up to the multiplicities described by the hard matrix element at NLO. The two b-jet bin in the 5FS -described only at LO by the hard matrix element -is affected by larger discrepancies between the two Monte Carlos, which exceed the uncertainty bands. Larger b-jet multiplicities, generated only at the shower level via gluon splittings, are affected by discrepancies as large as 100%: Herwig++ predicts less b jets than Pythia8.
Generally speaking, the difference between the Monte Carlos is smaller for the 4FS tha n for 5FS predictions. Indeed, the 4FS has more differential information at the matrixelement level, which reduces the effects of the shower. The only case worth mentioning, in which the Monte Carlo dependence is larger in the 4FS than in the 5FS is the distance between two hardest b jets at large separations. This observable is closely related to the distance of the two hardest B hadrons when two b jets are required (see right panel in figure 19). For both observables, at large separations, the 4FS predictions matched with Herwig++ lie very close to the 5FS NLO+PS predictions, while the matching with Pythia8 yields a higher tail. However, we point out that such discrepancy between Her-wig++ and Pythia8 in the 4FS barely exceeds the scale uncertainty bands.
Finally, we remark that similar conclusions can be drawn for the heavier charged Higgs case for most of the studied observables. However, a heavier charged Higgs reduces the differences between the two Monte Carlos (in particular in the 5FS) in the transverse momentum distribution of the two hardest B hadrons/b jets.

Conclusions
We have presented predictions for the production of a heavy charged Higgs boson in a type-II 2HDM, by explicitly considering a lighter Higgs scenario (m H − = 200 GeV) and a heavier one (m H − = 600 GeV). The only parameters that enter the calculation are the particle masses and tan β, thus the results are rather generic. Our predictions have been presented for tan β = 8 in a type-II 2HDM, but, through simple rescaling of the Yukawas,   they are applicable to any tan β value, to type-I 2HDMs (see section 6 of ref. [25] for details) and to the inclusion of the dominant SUSY corrections, e.g. in the MSSM.
For the first time, a fully differential computation has been performed in the 4FS at fNLO and NLO+PS accuracy. We have exploited the automatised Mad- Yukawa coupling in the MS scheme, which has the advantage with respect to the default on-shell scheme of resumming large logarithms of m b /µ R .
Our results indicate that a reduced shower scale with respect to the default one in MadGraph5 aMC@NLO improves the matching between parton shower and fixed-order results at large transverse momenta. For this scale choice, we discussed the effects of incorporating NLO corrections and the matching to the parton shower in the 4FS simulations by considering a number of differential observables. We found NLO corrections to be generally flat with our choice of the renormalisation and factorisation scales. Effects due to the parton shower are important in the Sudakov-dominated regions (jets and B hadrons at low p T ), or when observables sensitive to the b-quark fragmentation are considered. These observations have been made separately for the y 2 b and y 2 t terms. On the other hand, we argued that the y b y t contribution, appearing only in the 4FS, can be safely neglected, since its size is smaller than ∼ 5% of the total cross section (for tan β = 8 and m H − = 200 GeV) and is well within the scale uncertainty of the computation. For larger Higgs masses or different values of tan β it is further suppressed.
Besides discussing the new results in the 4FS, we provided a comprehensive comparison to the ones in the 5FS, consistently generated within MadGraph5 aMC@NLO. The inclusion of NLO(+PS) corrections in the two schemes improves their mutual agreement at the level of shapes. This agreement follows, although to a minor extent, from the reduced shower scale choice. Differences remain, however, and they are particularly sizeable for observables related to b jets and B hadrons. Given these differences, it was vital to carry out an unbiased analysis of our results in order to acquire the most reliable predictions for this class of observables. The proper simulation of the signal will be crucial for the experiments to fully exploit the potential of the data collected in charged Higgs searches at the LHC.
Our final recommendation is to use 4FS predictions for any realistic signal simulation in experimental searches. This recommendation is backed by two sets of evidences: first we have proven that, for a large number of observables, the 4FS prediction provides a better description of the final state kinematics; second, it reduces the systematic error related to the usage of a given parton shower. Moreover, when matching the NLO calculation to the shower, we recommend to use a lower shower scale (by setting F = 4 in our case). This choice provides a better matching to the fixed-order computation at large transverse momenta, slightly reduces the parton shower dependence and improves the agreement of four-and five-flavour scheme computations.
Any user interested in the simulation of charged Higgs production with Mad-Graph5 aMC@NLO is strongly encouraged to contact the authors. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.