Associated production of charged Higgs bosons and top quarks with POWHEG

The associated production of charged Higgs bosons and top quarks at hadron colliders is an important discovery channel to establish the existence of a non-minimal Higgs sector. Here, we present details of a next-to-leading order (NLO) calculation of this process using the Catani-Seymour dipole formalism and describe its implementation in POWHEG, which allows to match NLO calculations to parton showers. Numerical predictions are presented using the PYTHIA parton shower and are compared to those obtained previously at fixed order, to a leading order calculation matched to the PYTHIA parton shower, and to a different NLO calculation matched to the HERWIG parton shower with MC@NLO. We also present numerical predictions and theoretical uncertainties for various Two Higgs Doublet Models at the Tevatron and LHC.


I. INTRODUCTION
One of the most important current goals in high-energy physics is the discovery of the mechanism of electroweak symmetry breaking. While this can be achieved, as in the Standard Model (SM), with a single Higgs doublet field, giving rise to only one physical neutral Higgs boson, more complex Higgs sectors are very well possible and in some scenarios even necessary. E.g., in the Minimal Supersymmetric SM, which represents one of the most promising theories to explain the large hierarchy between the electroweak and gravitational scales, a second complex Higgs doublet is required by supersymmetry with the consequence that also charged Higgs bosons should exist.
At hadron colliders, the production mechanism of a charged Higgs boson depends strongly on its mass. If it is sufficiently light, it will be dominantly produced in decays of top quarks, which are themselves copiously pair produced via the strong interaction. Experimental searches in this channel have been performed at the Tevatron by both the CDF [1] and D0 [2] collaborations and have led to limits on the top-quark branching fraction and charged Higgs-boson mass as a function of tan β, the ratio of the two Higgs vacuum expectation values (VEVs), for various Two Higgs Doublet Models (2HDMs). However, if the charged Higgs boson is heavier than the top quark, it is dominantly produced in association with top quarks with a semi-weak production cross section. The D0 collaboration have searched for charged Higgs bosons decaying into top and bottom quarks in the mass range from 180 to 300 GeV and found no candidates [3]. At the LHC, the ATLAS (and CMS) collaborations have already excluded top-quark branching ratios to charged Higgs bosons with masses of 90 (80) to 160 (140) GeV and bottom quarks above 0.03−0.10 (0.25−0.28) using 1.03 fb −1 (36 pb −1 ) of data taken at √ S = 7 TeV [4,5]. At 14 TeV and with an integrated luminosity of 30 fb −1 , the discovery reach may be extended to masses of about 600 GeV using also the decay into tau leptons and neutrinos [6,7]. It may then also become possible to determine the spin and couplings of the charged Higgs boson, thereby identifying the type of the 2HDM realized other new physics sources [10].
In this paper, we concentrate on the associated production of top quarks and charged Higgs bosons at hadron colliders, which is of particular phenomenological importance for a wide range of masses and models. Conversely, s-channel single, pair, and associated production of charged Higgs bosons with W -bosons are less favorable in most models. Isolation of this signal within large SM backgrounds, e.g. from top-quark pair and W -boson associated production, and an accurate determination of the model parameters require precise predictions that go beyond the next-to-leading order (NLO) accuracy in perturbative QCD obtained previously [11][12][13]. We therefore present here details of our re-calculation of this process at NLO using the Catani-Seymour (CS) dipole subtraction formalism [14,15]. The virtual loop and unsubtracted real emission corrections are then matched with a parton shower (PS) valid to all orders in the soft-collinear region using the POWHEG method [16,17] in the POWHEG BOX framework [18]. A similar calculation has been presented using the Frixione-Kunszt-Signer (FKS) subtraction formalism [19] and matching to the HERWIG PS with the MC@NLO method [20]. Other new physics processes recently implemented in MC@NLO include, e.g., the hadroproduction of additional neutral gauge bosons [21]. Unlike MC@NLO, POWHEG produces events with positive weight, which is important when the experimental analysis is performed via trained multivariate techniques. POWHEG can be easily interfaced to both HERWIG [22] and PYTHIA [23] and thus does not depend on the Monte Carlo (MC) program used for subsequent showering.
The remainder of this paper is organized as follows: In Sec. II, we present details of our NLO calculation of top-quark and charged Higgs-boson production. We emphasize the renormalization of wave functions, masses, and couplings, in particular the one of the bottom Yukawa coupling, in the virtual loop amplitudes as well as the isolation and cancellation of soft and collinear divergences with the Catani-Seymour dipole formalism in the real emission amplitudes. The implementation in POWHEG is described in Sec. III. As the associated production of top quarks with charged Higgs bosons is very similar to the one with W -bosons [24], we concentrate here on the differences of the two channels. We also emphasize the non-trivial separation of the associated production from top-quark pair production with subsequent top-quark decay in scenarios, where the charged Higgs boson is lighter than the top quark, using three methods: removing completely doubly-resonant dia- FIG. 1: Tree-level diagrams for the associated production of charged Higgs bosons and top quarks at hadron colliders in the s-channel S and the t-channel T .
grams, subtracting them locally in phase space, or including everything, so that top-quark pair production with the subsequent decay of an on-shell top quark into a charged Higgs boson is effectively included at leading order (LO). In Sec. IV, we present a detailed numerical comparison of the new POWHEG implementation to the pure NLO calculation without PS [12], to a tree-level calculation matched to the PYTHIA parton shower [25], and to the MC@NLO implementation with the HERWIG PS [20]. We also give numerical predictions and theoretical uncertainties for various 2HDMs at the Tevatron and LHC. Our conclusions are summarized in Sec. V.

A. Organization of the calculation
At the tree level and in the five-flavor scheme with active bottom (b) quarks as well as gluons (g) in protons and antiprotons, the production of charged Higgs bosons (H − ) in association with top quarks (t) occurs at hadron colliders via the process b(p 1 ) + g(p 2 ) → H − (k 1 ) + t(k 2 ) through the s-and t-channel diagrams shown in Fig. 1 A NLO calculation in a four-flavor scheme, where the bottom quark is treated as massive and generated by the splitting of an initial gluon, has been presented elsewhere [26], but the effect of the bottom mass through the parton densities was subsequently found to be strongly suppressed compared to its impact on the bottom Yukawa coupling [27]. The fourmomenta of the participating particles have been ordered in accordance with the POWHEG scheme, where the initial-state particles with four-momenta p 1 and p 2 are followed by the four-momentum k 1 of the final-state massive colorless particle and then the four-momentum k 2 of the outgoing massive colored particle. The additional radiation of a massless particle (gluon or light quark), that occurs at NLO in real emission diagrams, is assigned the last four-momentum k 3 .
The hadronic cross section is obtained as usual as a convolution of the parton density functions (PDFs) f a/A, b/B (x a,b , µ 2 F ) with the partonic cross section with partonic center-of-mass energy s = x a x b S, S being the hadronic center-of-mass energy.

Its LO contribution
is obtained from the spin-and color-averaged squared Born matrix elements |M Born | 2 through integration over the two-particle phase space dΦ (2) and flux normalization.

B. Virtual corrections and renormalization scheme
Like the LO partonic cross section σ LO , its NLO correction has a two-body final-state contribution which consists of the virtual cross section dσ V , i.e. the spin-and color-averaged interference of the Born diagrams with their one-loop corrections, and the Born cross section dσ LO convolved with a subtraction term I, which can be written as 2 H, t; b, g | I(ǫ) | H, t; b, g 2 and which removes the infrared singularities present in the virtual corrections. The three-body final-state contribution σ N LO{3} and the finite remainders σ N LO{2} (x, ...) of the initial-state singular terms will be described in the third part of this section.
The ultraviolet divergencies contained in the virtual cross section dσ V have been made explicit using dimensional regularization with D = 4 − 2ǫ dimensions and are canceled against counterterms originating from multiplicative renormalization of the parameters in the Lagrangian. In particular, the wave functions for the external gluons, bottom and top quarks are renormalized in the MS scheme with where ∆ U V = 1/ǫ − γ E + ln 4π, γ E is the Euler constant, N C = 3 and N F = 6 are the total numbers of colors and quark flavors, respectively, and C F = (N 2 C − 1)/(2N C ). The counterterm for the strong coupling constant α S = g 2 S /(4π) is computed in the MS scheme using massless quarks, but decoupling explicitly the heavy top quark with mass m t from the running of α S [28]. The top-quark mass entering in the kinematics and propagators is renormalized in the on-shell scheme, On the other hand, we perform the renormalization of both the bottom and top Yukawa couplings in the MS scheme, This enables us to factorize the charged Higgs-boson coupling at LO and NLO, making the QCD correction (K) factors independent of the 2HDM and value of tan β under study. In particular, in Eq. (13) we do not subtract the mass logarithm, but rather resum it using the running quark massesm (14) in the Yukawa couplings, where for m b < µ R < m t and for µ R > m b,t [29]. The starting values of the MS masses are obtained from the on-shell masses M Q through the relation with K b ≈ 12.4 and K t ≈ 10.9 [30,31].
After the renormalization of the ultraviolet singularities has been performed as described above, the virtual cross section contains only infrared poles. These are removed with the second term in Eq. (8), i.e. by convolving the Born cross section with the subtraction term [14,15] where in our case I 2 (ǫ, µ 2 ; {k 2 , m t }) = 0, since there are no QCD dipoles with a final state emitter and a final state spectator. The dipoles depending on one initial-state parton (a = b, g) with four-momentum p i (i = 1, 2) are where T a,t denotes the color matrix associated to the emission of a gluon from the parton a or the top quark t, the dimensional regularization scale µ is identified with the renormalization scale µ R , and s ta = s at = 2p i k 2 . The kernels consist of the singular terms with Q 2 ta = Q 2 at = s ta + m 2 t + m 2 a and the non-singular terms The constant κ in Eq. (19) is a free parameter, which distributes non-singular contributions between the different terms in Eq. (7). The choice κ = 2/3 considerably simplifies the gluon kernel. For massive quarks, one has in addition while and with T R = 1/2 and N f = 5 the number of light quark flavors. The last term in Eq. (18) depends on both initial-state partons. Since the process we are interested in involves only three colored particles at the Born level, the color algebra can be performed in closed form.
To be concrete, we have

C. Real corrections
The second term in Eq. (7) includes the spin-and color-averaged squared real emission matrix elements |M 3,ij (k 1 , k 2 , k 3 ; p 1 , p 2 )| 2 with three-particle final states and the corresponding unintegrated QCD dipoles D, which compensate the integrated dipoles I in the previous section.
Both terms are integrated numerically over the three-particle differential phase space dΦ (3) .
The real emission processes can be grouped into the four classes where the second process (b) can be obtained from the first one (a) by crossing the fourmomenta k 3 and −p 1 and multiplying the squared matrix element by a factor of (−1) to take into account the crossing of a fermion line. The processes in the two other classes (c) and (d) can interfere when q = b, but these contributions are numerically negligible due to the comparatively small bottom quark parton distribution function. Process (d) is furthermore convergent for q = u, d, s and c.
The sum over the dipoles in Eq. (35) includes initial-state emitters ab with both initialand final-state spectators c (D ab,c and D ab c ) and the final-state emitter ab with initial-state spectators c (D c ab ). For the three divergent processes, we have (a) : Denoting by a the original parton before emission, b the spectator, and i the emitted particle, the dipole for initial-state emitters and initial-state spectators is given by where the momentum of the intermediate initial- The necessary splitting functions V ai,b for {ai, b} = {qg, g; gg, q; gq, g; qq, q} can be found in Ref. [14]. The dipole for initial-state emitters and a final-state spectator, which is in our case the top quark t, is given by where the momentum of the intermediate initial- The necessary splitting functions V ai t for {ai, t} = {qg, t; gg, t; gq, t; qq, t} can be found in Ref. [15]. Finally, the dipole for final-state emitter (the top quark t) and initial-state spectator a is given by where the momentum of the initial parton a is shifted top µ The required splitting function V a gt can again be found in Ref. [15]. The last terms in Eq. (7) are finite remainders from the cancellation of the ǫ-poles of the initial-state collinear counterterms. Their general expressions read and similarly for (a ↔ b) and (p 1 ↔ p 2 ). The color-charge operators K and P are explicitly given in Ref. [15].

III. POWHEG IMPLEMENTATION
The calculation in the previous section has been performed using the Catani-Seymour dipole formalism for massive partons [14,15]. For the implementation of our NLO calculation in the POWHEG Monte Carlo program, we need to retain only the Born process, the finite terms of the virtual contributions, and the real emission parts of our calculation, since all necessary soft and collinear counterterms and finite remnants are calculated automatically by the POWHEG BOX in the FKS scheme [19]. Soft and collinear radiation is then added to all orders using the Sudakov form factor. In this section, we briefly describe the three relevant contributions, following closely the presentation in Ref. [18], and address the nontrivial separation of the associated production of charged Higgs bosons and top quarks from top-quark pair production with subsequent top-quark decay in scenarios, where the charged Higgs boson is lighter than the top quark.

A. Born process
In the POWHEG formalism, a process is defined by its particle content. Each particle is encoded via the Particle Data Group numbering scheme [32] except for gluons, which are assigned the value zero. The order of the final state particles has to be respected. Colorless particles are listed first, then heavy colored particles, and finally massless colored particles.
The Born processes (two with respect to the different bg and gb initial states) are defined with flst nborn = 2 and are listed as and in the subroutine init processes.
In the subroutine born phsp, the integration variables xborn(i) for the Born phase space are generated between zero and one. The hadronic cross section is then obtained from the differential partonic cross section dσ via the integration (see Eq. (4)) where f i/I is the PDF of parton i inside hadron I with momentum fraction x i and where we have performed the change of variables The integration limits are given in Tab. I. The Jacobian for the change of integration variables from xborn(i) to (τ, y, t) has to be multiplied with 2π for the integration over the azimuthal angle φ, which is randomly generated by POWHEG. The different kinematical variables can then be constructed in the where k is to be varied around two for uncertainty studies. Both the born phsp and the set fac ren scales subroutines can be found in the file Born phsp.f.
before summing over the initial gluon polarizations as well as where g µν is the metric tensor.

B. Virtual loop corrections
The renormalized virtual cross section is defined in dimensional regularization and in the POWHEG convention by where |M Born | 2 is now the squared Born matrix element computed in D = 4 − 2ǫ dimensions and where the remaining double and simple infrared poles are proportional to The POWHEG implementation needs then only the finite coefficient V f in , which has been At this point the problem arises how to separate the two production mechanisms. In the literature two methods have been proposed: Diagram Removal (DR) and Diagram Subtraction (DS) [33]. Both remove the top-quark resonance from the cross section, but the procedure for combining top pair production with the associated production is not completely clear. If we separate the amplitudes of a real emission process with colliding partons a and b into contributions M tt ab , which proceed through tt-production, and contributions M tH − ab , which do not, squaring the amplitudes gives rise to three different quantities: The term D ab contains neither collinear nor soft singularities, while the interference term I ab contains integrable infrared singularities. These terms are therefore sometimes referred to as subleading with respect to those in S ab , which contains all infrared singularities and must be regularized, e.g., via the subtraction formalism. DR requires removing tt production at the amplitude level. The only contributing element is then S ab . Since it contains all divergencies, the dipoles used in the m H > m t case remain valid. In the DS scheme, one subtracts from the cross section the quantity locally in phase space. The momenta are reorganized so as to put thet quark on its mass shell. Although gauge invariant, this procedure is still somewhat arbitrary. We therefore introduce here a third option, where nothing is removed or subtracted from the associated production, but simply the full production cross section is retained. Once a sample of events is generated, one can then still decide to remove events near the resonance region and replace them with events obtained, for example, with a full NLO implementation of tt production.
In our POWHEG code, we implemented all three methods described above. DR is the simplest case. If the flag DR is set to one in the file powheg.input, the resonant diagrams of While with the DR or the full scheme the fraction of negative weights is very small, this is not the case in the DS scheme. Here the real cross section can become negative in certain kinematic regions, so that POWHEG must then be run with the flag withnegweights set to one. Negatively weighted events are then kept, but are hard to interpret, since they correspond to the subtraction of an ad hoc quantity from the cross section.
Removing diagrams at the amplitude level causes the loss of gauge invariance. A considerable part of Ref. [33] has been dedicated to the analysis of the corresponding impact on W t production. There, different gauges were considered for the gluon propagator, and differences at the per-mille level were found. Note, however, that gauge invariance is not only spoiled through the gluon propagator, but also when the polarization sum of external gluons is replaced by for simplicity. Here, k µ is the four-momentum, λ is the polarization, and ǫ µ (k, λ) is the polarization vector of the external gluon. Eq. (59) includes not only physical transverse, but also non-physical gluon polarizations that must be canceled by ghost contributions. Removing individual diagrams then causes the loss of gauge invariance. We therefore abandon the use of the simple polarization sum, Eq. (59), and sum instead only over physical states with where η µ is an arbitrary four-vector transverse to the polarization vector ǫ µ . When calculating a gauge invariant quantity, the η-dependence would drop out, but this will not be the case in DR as argued above. For the channels with two external gluons and incoming four-momenta p 1 and p 2 , we choose for the polarization vectors  [37]. Both hypotheses have recently been compared with the result that the latter appears to be more stable under quantum corrections, but that the two hypotheses are largely equivalent at tree level [38].
In a general 2HDM, one introduces two complex SU(2)-doublet scalar fields where the Vacuum Expectation Values (VEVs) v 1,2 of the two doublets are constrained by the W -boson mass through v 2 = v 2 1 +v 2 2 = 4m 2 W /g 2 = (246 GeV) 2 [39]. The physical charged Higgs bosons are superpositions of the charged degrees of freedom of the two doublets, In the Type-I 2HDM, only Φ 2 couples to the fermions in exactly the same way as in the minimal Higgs model, while Φ 1 couples to the weak gauge bosons [40]. The Feynman rules where V ij is the CKM matrix and P L,R = (1∓γ 5 )/ √ 2 project out left-and right handed quark eigenstates. As can be seen from Tab. III, these couplings are the same in the lepton-specific 2HDM.
In the Type-II 2HDM, Φ 2 couples to up-type quarks and Φ 1 to down-type quarks and charged leptons. The Feynman rules for charged Higgs-boson couplings to quarks in this model, with all particles incoming, are As can again be seen from Tab. III, they are identical to those in the flipped 2HDM [41].
Since the NFC and MFV hypotheses allow for the possibility that the two Higgs doublets couple to quarks with arbitrary coefficients A i u,d , there exists also the possibility of a more general 2HDM, sometimes called Type-III 2HDM [42]. In this case, the Feynman rules for the charged Higgs-boson couplings to quarks, with all particles incoming, are where the family-dependent couplings A i u,d read  [42]. Since general colorsinglet Higgs-boson couplings and (theoretically possible) color-octet Higgs bosons induce different QCD corrections, we will not study these scenarios numerically. In the literature, one may also find Type-III 2HDMs where no flavor symmetry is imposed and FCNCs are avoided by other methods, e.g. by the small mass of first and second generation quarks [43].
These models then allow for the couplings of charged Higgs bosons to bottom and charm quarks, which induces a phenomenology that is different from the one studied in this paper.

C. Predictions for various 2HDMs
The calculation presented in the previous sections was performed in a generic way, which makes it possible to use the result for various models with charged Higgs bosons. Out of the models mentioned in the last section, our calculation is in particular valid for the Type-I and Type-II 2HDMs. In this subsection, we perform a numerical analysis for a set of typical collider scenarios for these 2HDMs. As in these models the scattering matrix element is directly proportional to the Higgs-top-bottom quark coupling even at NLO, the type of the model has no influence on kinematic distributions apart from their normalization to the total cross section.
We therefore concentrate here on the total cross sections and on the uncertainties both In all scenarios, both at the Tevatron and at the LHC, the next-to-leading order correction is substantial, ranging from 57% at the Tevatron in the Type-I 2HDM to 38% at the LHC in the same model. Apart from enhancing the total cross section, including the NLO correction At leading order, the strong scale dependence comes from the strong coupling constant and from the Yukawa coupling in the tree-level amplitude. Including higher-order corrections, this uncertainty is dramatically reduced in some scenarios.
Another large source of error stems from the parton distribution functions. We use the CT10 NLO PDF set with its error PDF sets to determine the error coming from the uncertainty contained in determining the parton content of the colliding hadrons. The process considered here is extremely sensitive to the gluon distribution function through having a gluon in the initial state and through having a heavy-quark initial state, which is radiatively generated from the gluon PDF. Moreover, the production of a heavy Higgs boson in association with a top quark probes the higher x content of the initial-state (anti-)proton.
The values of Bjorken-x probed can be expressed as which at the Tevatron leads to typical values of x ∼ 0.3. This is exactly the region where the gluon PDF is poorly known, which translates into large PDF uncertainties on the cross section at the Tevatron. At the LHC, the Bjorken-x probed is x ∼ 0.1, and the PDF uncertainties are therefore much smaller.

D. Checks of the NLO calculation and comparisons with POWHEG
As a check of the numerical implementation of our analytical results, we have compared our complete NLO calculation obtained with the Catani-Seymour dipole formalism with the one performed previously with a phase-space slicing method using a single invariantmass cutoff [12], which had in turn been found to agree with a calculation using a two (soft and collinear) cutoff phase-space slicing method [11]. We found good agreement for all differential and total cross sections studied, but refrain from showing the corresponding figures here, since the fixed-order results are well-known.
For the remainder of the analysis, we will constrain ourselves to the Type-II 2HDM, as the kinematic distributions have the same features in both Type-I and Type-II 2HDMs. In all of our discussion, we consider three collider scenarios: • Tevatron, √ S = 1.96 TeV, • LHC, √ S = 7 TeV, and   This provides a good consistency check of the matching procedure. Before schemes to match parton showers with full NLO calculations were developed, the importance of the contribution of this particular two-to-three process and the perturbative origin of the b-quark density had already been recognized [25]. It had been proposed to supplement the LO calculation by this particular two-to-three process and to remove the overlap by subtracting the doubly counted (DC) term where with P qg (z) the g → q splitting function, f g (x, µ 2 F ) the gluon PDF, and z the longitudinal gluon momentum fraction taken by the b-quark. The two-to-three and double-counting processes had been implemented in an addition to PYTHIA called MATCHIG.
With our full NLO calculation matched to PYTHIA within the POWHEG BOX, it is now possible to compare the two approaches numerically. The results are shown in Fig. 10. Since the normalization of the MATCHIG prediction is still effectively of LO, we have normalized all distributions to their respective total cross sections in order to emphasize the shapes of the distributions. One observes that when both the MATCHIG (black) and POWHEG (red) predictions are matched to the PYTHIA parton shower, there is very little difference, even at low p T and large ∆φ of the top-Higgs pair. Only at large p T and small ∆φ the differences become sizable, which can be attributed to the fact that MATCHIG includes only one of the four classes of real-emission processes, while our POWHEG prediction includes also the quark-initiated real-emission processes. Let us emphasize again that while the spectra are already quite well described with MATCHIG, their normalization is only accurate to LO and not NLO as in POWHEG.

G. Comparison with MC@NLO
In a recent publication, two of us and a number of other authors have matched a NLO calculation performed with the FKS subtraction formalism to the HERWIG PS with the MC@NLO method [20]. It is therefore mandatory that we compare in this paper this previous work with our new POWHEG implementation, which we do in Fig. 11. Note that here we employ a value of tan β = 30 as in the MC@NLO publication. In both calculations, we use the HERWIG PS in order to emphasize possible differences in the matching methods and not those in the parton shower. We also normalize the differential cross sections again to the total cross section for a better comparison of the shapes of the distributions.
As in the other comparisons, the rapidity distributions of the charged Higgs boson (top right) and the top quark (center right) show little variation, confirming the consistency of the two calculations. However, the corresponding p T -spectra (top and center left) are slightly harder with the MC@NLO matching than in POWHEG. This behaviour is known from other processes [24,44]. It is less pronounced in the p T -distribution of the top-Higgs pair, shown on a logarithmic scale (bottom left). Since we are using the HERWIG PS, the rise at small azimuthal angle ∆φ (bottom right) is not very strong with MC@NLO and only slightly more so with POWHEG. In total, all of these differences are similarly small in the production of a top quark with a W -boson [24] and with a charged Higgs boson at the LHC.

H. Diagram Removal, Diagram Subtraction, and no subtraction
If the charged Higgs boson was lighter than the top quark, it would dominantly be created in top-pair production and the decay of an (anti-)top quark into it. As discussed above, one must then find a suitable definition to separate this process from the associated top-Higgs production discussed in this paper. In addition to the Diagram Removal (DR) and Diagram Subtraction (DS) methods discussed above, we introduce here also the option of not removing or subtracting anything from the associated production, but simply retaining the total production cross section, which then allows for the removal of fully simulated events near the resonance region and replacing them with events obtained, e.g., with a full NLO implementation of tt production. The results are shown in Fig. 12. The theoretical pros and cons and the numerical differences of Diagram Removal and Diagram Subtraction have been discussed extensively above and also elsewhere [20]. It is clear from Fig. 12 that the numerical difference of DR vs. DS is much less pronounced than the difference of both with respect to no removal or subtraction at all. We emphasize that the total cross section is continuous across the m H = m t threshold in all three schemes (see also Ref. [27]).
The differences of POWHEG and MC@NLO are small for m H < m t in both the DR and DS schemes, as can be seen when comparing Figs. 13 and 14. This coincides nicely with our observation above that these differences should be as small as in the associated production of W -bosons and top quarks [24].

V. CONCLUSION
In this paper, we presented a new NLO calculation of the associated production of charged Higgs bosons and top quarks at hadron colliders using the Catani-Seymour dipole subtraction formalism and matched it to parton showers with the POWHEG method. We discussed the It will now be very interesting to observe the impact of our work on the experimental search for charged Higgs bosons. The numerical code and technical support is, of course, available from the authors.