Higgs production in association with bottom quarks

We study the production of a Higgs boson in association with bottom quarks in hadronic collisions, and present phenomenological predictions relevant to the 13 TeV LHC. Our results are accurate to the next-to-leading order in QCD, and matched to parton showers through the MC@NLO method; thus, they are fully differential and based on unweighted events, which we shower by using both Herwig++ and Pythia8. We perform the computation in both the four-flavour and the five-flavour schemes, whose results we compare extensively at the level of exclusive observables. In the case of the Higgs transverse momentum, we also consider the analytically-resummed cross section up to the NNLO+NNLL accuracy. In addition, we analyse at ${\cal O}(\alpha_S^3)$ the effects of the interference between the $b\bar{b}H$ and gluon-fusion production modes.


Introduction
Data collected at the LHC so far fully support the hypothesis that the observed resonance with a mass of about 125 GeV is the scalar boson predicted by the Brout-Englert-Higgs symmetry breaking mechanism [1,2] of SU (2) L × U (1) Y , as implemented in the Standard Model (SM) [3][4][5]. In such a theory, the strengths of the Higgs boson couplings to all the elementary particles (including the Higgs itself) are universally set by the corresponding masses. A global fit to the production rates, that employs several different final states, shows a 10-20% agreement with SM predictions for the measured couplings to fermions and to vector bosons [6][7][8]. Conversely, no information has yet been obtained on the Higgs self-coupling.
To probe elementary couplings in both decay and production is not only interesting, but also quite useful in order to reduce the theoretical as well as the experimental uncertainties. The case of the Higgs coupling to bottom quarks (y b ) is rather special in this respect: while its strength is significantly smaller than those relevant to vector bosons and to top quarks, for a Higgs mass of about 125 GeV the H → bb mode dominates the total decay width, owing to phase-space factors. Unfortunately, this is not particular helpful from an experimental viewpoint, for several reasons: the total width is extremely small in absolute value; the backgrounds which feature b quarks are immense; and the relative partial decay widths are difficult to determine with some accuracy, precisely because the H → bb branching ratio is close to one. A Higgs production mode that features a bbH coupling is thus a more viable alternative. There are two main such modes: the dominant one is a contribution to gluon fusion, where bottom-quark loop amplitudes interfere with top-quark loop amplitudes; the second-largest is associated production (bbH henceforth) -with this, one understands processes that at the Born level receive contributions from tree graphs that include a b-quark line which radiates a Higgs boson. The b-quark contributions to gluon fusion are of the order of a few percent for the total rate, and up to 10% for a Higgs produced at small transverse momentum (see e.g. refs. [9][10][11][12]). Conversely, the fully-inclusive bbH cross section is much smaller: it is of the same order as the ttH one (around 0.5 pb), i.e. the fifth largest after gluon fusion, VBF, W H, and ZH. However, this inclusive rate decreases dramatically when conditions are imposed (which means minimal transverse momentum and centrality requirements, typical of b-tagging) that allow one to render it distinguishable from other production mechanisms.
The SM picture outlined above might be significantly modified by beyond-SM effects. For example, in extended Higgs sectors the couplings of the scalar particle(s) to bottom quarks may be enhanced, typically by a factor 1/ cos β or tan β for pseudoscalars in generic or supersymmetric 2HDM's. Given that a scalar sector richer than that of the SM has not yet been ruled out experimentally, this is a fact that one must bear in mind, and that constitutes a strong motivation for pursuing the study of scalar-particle production in association with b quarks.
In addition, and regardless of phenomenological motivations, bbH production is interesting in its own right theoretically, and has in fact generated much discussion in the past. The main reason for this is that, as for all mechanisms that feature b quarks at the hard-process level, there are two ways of performing the computation, which are usually called four-and five-flavour schemes (abbreviated with 4FS and 5FS henceforth, respectively). These two schemes are supposed to address issues that arise in different kinematic regimes, which one may classify by using a hard scale (that we denote by Q) typical of the process. Lest we complicate the discussion with a proliferation of scales, for this discussion we assume to be dealing with fully-inclusive observables.
Four-flavour scheme computations are relevant to those cases where the physical mass of the b quark (m b ∼ 5 GeV) plays the role of a hard scale, which means: (1.1) In the context of a factorisation theorem, this implies that the b quark must be treated as a massive object at the level of the short-distance cross section, where it never appears in the initial state; that the factorisation formula neglects terms of order Λ QCD /m b ; and that observables with tagged final-state b quarks (in association with whatever other object) can be computed in perturbation theory. In the case at hand, the leading order (LO) partonic processes are therefore: gg → bbH , qq → bbH , (1.2) with q a light quark. The condition m b ∼ Q is sufficient for eq. (1.1) to be satisfied, but it is not necessary. In general, by computing the cross section in perturbation theory, at any perturbative order the 4FS result will feature terms proportional to α k S log k (m b /Q). These are harmless when m b ∼ Q, but if Q ≫ m b they might spoil the "convergence" of the perturbative series, in spite of the condition in eq. (1.1) being still fulfilled. This is hardly an unusual situation, and certainly one which is not peculiar to b physics: when a series is not well-behaved, one re-organises its expansion, by resumming the appropriate logarithmic terms. On the other hand, the relative smallness of m b renders the presence of large logarithms a likely occurrence. If the crucial characteristic of an observable is that of being dominated by such logarithms, then eq. (1.1) is not really relevant, and a more appropriate description of the dominant kinematic regime is: which is what the five-flavour scheme computations deal with. The easiest way to achieve eq. (1. 3), and one that lends itself particularly well to perturbative computations, is that of setting m b = 0 at the level of short-distance cross sections. In a factorisation theorem, thus, the b quark will be treated on equal footing as all of the other light quarks including, in particular, the fact that it could appear in the initial state. This implies that, in the 5FS the LO partonic process for Higgs production in association with b quarks is: The fact that m b = 0 has to be regarded as a technical mean to achieve the resummation of large logarithms, which resummation is the reason why the 5FS has been introduced in the first place. In the case of a fully-inclusive cross section, the logarithms are indeed effectively resummed through the Altarelli-Parisi evolution equations relevant to the b-quark PDF 1 . From the purely theoretical viewpoint, 5FS computations have the advantage of being much simpler than their 4FS counterparts: the process in eq. (1.4) is a 2 → 1 O(α 0 S ) one, while those in eq. (1.2) are 2 → 3 and of O(α 2 S ). This advantage comes at a steep price: the calculation associated with eq. (1.4) does provide one with a much more limited information than eq. (1.2), owing to the absence of final-state b quarks in the former. In the 5FS, such information can only be recovered by considering higher orders in α S , and in so doing the 5FS starts losing the advantage mentioned above. Explicitly for the case of bbH production, while in the 4FS an O(α 2 S ) term is always perturbatively leading, regardless of whether one considers 0-, 1-, or 2-b tag observables, in the 5FS the leading contribution is O(α S ) for 1-b tag observables, and O(α 2 S ) for 2-b tag observables. Furthermore, in the context of a 5FS computation a b-tagged object leads to unphysical results if defined with too small a p T (unphysical since the b-tagged cross section is larger than the fully-inclusive one; it diverges at p T → 0), owing to the mass of the b quark having been set equal to zero; this problem does not affect 4FS results. Finally, b-tagged objects in the 5FS cannot coincide with b quarks (which, conversely, is the case in the 4FS), because the corresponding cross section is not finite order-by-order in perturbation theory: the b quarks must either be integrated over, or be part of jets, or be converted into b-flavoured hadrons through the convolution with fragmentation functions.
One must also bear in mind that, in general, when considering differential observables new mass scales become relevant to the problem, which rapidly render the generalisation of eqs. (1.1) and (1.3) impractical, and make it difficult to decide a priori which calculational scheme is best suited to tackle the problem at hand.
While 4FS results lack logarithmic terms beyond the first few, 5FS results lack powersuppressed terms (m b /Q) n . Which of the two classes of terms is more important depends on the observable studied, that determines the dominant kinematic regime. In order to avoid the problems that this fact entails, schemes [13][14][15][16][17][18] have been proposed that allow one to combine, possibly in a systematic manner in perturbation theory, the logarithmic and power-suppressed terms in a single cross section, which is appropriate to all kinematic regimes. It is customary (but not necessary) to view these approaches as power-termsimproved 5FS calculations, which is sort of natural when one looks at fairly inclusive observables, for which one expects logarithms to be more important than power-suppressed terms, so that the latter are small corrections to the former. In the following and for the sake of clarity, by 5FS results we shall understand those that, at the level of short-distance cross sections, do not contain any power-suppressed terms.
Given what has been said so far, a typical rationale is the following: if logarithms are large, the 5FS should be superior to the 4FS; if they are not, and thus power-suppressed terms might be important, then 4FS approaches should be preferred. Unfortunately, the relative weight of logarithmic and power-suppressed terms is observable dependent (which is also complicated by the fact that a given observable can be potentially associated with different powers of α S in the four-and five-flavour schemes, as discussed above). However, one expects that, for processes and in regions of the phase space where both resummation and mass effects are not dominant, the two approaches should give similar results. As a matter of fact, at least for inclusive quantities 5FS and 4FS results are indeed generally similar (in particular, for bbH production), because of the two following facts. i) Logarithms are dominant, but they are especially so only at large Bjorken x's, and are always associated with phase-space suppression factors [19]. ii) At the lowest perturbative orders (LO and NLO), a reasonable agreement is found only by judiciously choosing the hard scales; in particular, arguments based on collinear dominance suggest that the optimal values of these scales in bbH production are significantly smaller than m H [20][21][22], i.e. than the hardness one would naively associate with this production process. As a consequence of these two facts, the few logarithms that appear in a 4FS fixed-order cross section approximate well numerically the leading logarithmic tower(s) present in the corresponding 5FS cross section (especially at the NLO), while power-suppressed terms are unimportant. There is an ample heuristic evidence of the facts above. However, one must keep in mind that such an evidence, and the arguments that support it, are chiefly relevant to inclusive variables. In a more exclusive scenario, it is important to assess the situation in an unbiased manner: this is one of the main aims of this paper.
What has been said so far has tacitly assumed fixed-order parton-level computations; when matrix elements are matched to parton showers (PS), some aspects of the picture do change. Let us start by considering 5FS approaches. For these, the foremost consequence of PS matching is the fact that even the O(α 0 S ) cross section of eq. (1.4), thanks to the backward evolution of the initial-state b's, will generate b-flavoured hadrons in the final state, which will render realistic any b-tagging requirements, regardless of the p T values where they are imposed. While this goes a long way towards improving 5FS predictions, one must not forget two facts. Firstly, it is well known that the backward evolution of massless b quarks is not trivial for Monte Carlos (MCs; see e.g. sect. 3.3 of ref. [23], and later here in sect. 3.3); this can lead to unexpected results for certain classes of observables, and to a marked MC dependence. In particular, the necessary kinematic reshuffling of the massless b quarks into massive b quarks can have a significant impact on the kinematics of final-state B mesons. Secondly, beyond LL the Altarelli-Parisi equations and the MC backward evolution do differ (and especially so for the treatment of b-quark thresholds); although generally small, these differences might become relevant in certain phase-space corners, where comparisons to data will help decide which description is best. As far as 4FS predictions are concerned, they are also improved by their matching with PS. In particular, small-p T initial-state emissions are Sudakov-suppressed, and effectively resummed by the parton showers 2 . Note, finally, that the matching to PS introduces in both the 4FS and 5FS extra power-suppressed effects, due to long-distance phenomena.
In summary, there are a number of interesting physics questions that are relevant to bbH production, which become especially crucial when exclusive quantities are studied, in particular in the context of computations matched to parton showers. The main results of this paper are the following.
• We present the first NLO computations matched to parton showers (NLO+PS) in the four-and five-flavour schemes. At the level of fixed-order NLO results (fNLO) we present fully-differential results which extend the scope of those available in the literature [24,25] in a very significant manner.
• We study, for the first time, the effect of the O(y b y t α 3 S ) interference term on differential distributions in the 4FS, in particular at the NLO+PS accuracy.
• We compare 4FS and 5FS NLO+PS predictions at the level of differential distributions, in order to further the arguments given above for such observables. For the inclusive Higgs transverse momentum, we also compare to the (N)NLO+(N)NLL analytical results of ref. [26].
All of our computations, bar those that feature analytical resummations, are performed in the automated MadGraph5 aMC@NLO framework [27].
The paper is organised as follows: in sect. 2 we report some generalities relevant to our 4FS and 5FS computations; phenomenological results are presented in sect. 3 -see in particular sect. 3.2 for 4FS predictions, and sect. 3.3 for 4FS-vs-5FS comparisons; we conclude in sect. 4.

Outline of the calculations
In this section, we briefly describe the physics contents relevant to our 4FS and 5FS predictions which have been calculated with MadGraph5 aMC@NLO. We remind the reader that MadGraph5 aMC@NLO contains all ingredients relevant to the computations of LO and NLO cross sections, with or without matching to parton showers. NLO results not matched to parton showers (called fNLO [27]) are obtained by adopting the FKS method [28,29] for the subtraction of the singularities of the real-emission matrix elements (automated in the module MadFKS [30]), and the OPP integral-reduction procedure [31] for the computation of the one-loop matrix elements (automated in the module MadLoop [32], which makes use of CutTools [33] and of an in-house implementation of the optimisations proposed in ref. [34] (OpenLoops)). Matching with parton showers is achieved by means of the MC@NLO formalism [35]. MadGraph5 aMC@NLO is maximally automated, since the only operations required by the user are to enter the process he/she is interested in generating, and the associated input parameters. The case of bbH production has one peculiarity related to the central role of the bottom Yukawa, which could not be handled by the public version available when ref. [27] was released, and which has been the object of a special treatment for the sake of this paper; more details are given in the next section.
Before introducing the features of our own computations, we briefly summarise the status of the bbH-production results available in the literature. As far as the 4FS is concerned, NLO fixed-order parton-level predictions have been presented in refs. [24,25], and later updated to the case of MSSM-type couplings [36], and to SUSY-QCD corrections in the MSSM [37]. The focus of these papers is the total cross section; only a handful of differential predictions have been shown. The literature is considerably richer for the 5FS, owing to its being technically simpler. NLO and NNLO QCD corrections for total rates were first computed in refs. [38,39] and in ref. [22], respectively. Differential parton-level predictions have been later made available: at the NLO for H + b and H+jet production [40,41], and at the NNLO for jet rates [42] and fully differential distributions [43]. The Higgs transverse momentum has been studied analytically at the NNLO in ref. [44], while resummed NLO+NLL and NNLO+NNLL results have been presented in ref. [45] and ref. [26], respectively. Prior to this paper, no NLO+PS predictions have been obtained in either scheme.

Four-flavour scheme
At the LO, the 4FS cross section receives contributions from the O(α 2 S ) 2 → 3 partonic subprocesses given in eq. (1.2); representative Feynman diagrams are displayed in fig. 1. The Higgs is always radiated off a b quark, and therefore the cross section at this perturbative order is proportional to y 2 b α 2 S . The coupling structure becomes more involved when one considers NLO corrections. As is well known, these can be classified as being of either virtual or real-emission origin; sample diagrams for these two classes are displayed in fig. 2 and fig. 3 respectively. Consider, in particular, the virtual diagrams of fig. 2a, 2b, and 2e: when the heavy quark that circulates in the loop is a top, the corresponding amplitude is proportional to y t , and does not feature the bottom Yukawa y b . All of the other diagrams, as well as those relevant to real emissions, have amplitudes proportional to y b . At the NLO, the cross section receives contributions from the interference of the one-loop diagrams with the Born ones, and from the squares of the real-emission diagrams. The squares of the one-loop diagrams, in turn, will enter the NNLO result. One can thus write the bbH 4FS cross section up to O(α 4 S ) as follows: (2.1) The bbH NLO results presented in the literature focus on the σ y 2 b term of eq. (2.1). We are not aware of existing predictions for σ y b yt at the level of differential observables, whose impact we shall discuss in sect. 3.2. Finally, all terms of O(α 4 S ) have been ignored here; note that at least those proportional to y 2 t are usually seen as NNLO contributions to the gluon-fusion cross section.
The fully automated MadGraph5 aMC@NLO [27] program can handle 4FS bbH production rather straightforwardly -the calculation is of a complexity similar to that of Zbb production, which could be studied [46] even with a version of the code much less powerful than the present one. However, the default treatment of Yukawa couplings in the code is that of an on-shell scheme renormalisation, which is not optimal in the case of bbH production, where the MS scheme has to be preferred [47]. Such a scheme is indeed what has been employed in previous NLO 4FS computations [24,25], since the use of an MS renormalized Yukawa y b (µ R ) has the advantage of resumming to all orders potentially large logarithms of m H /m b , when µ R ∼ m H is chosen. A change in the renormalisation scheme, and the UV counterterms it entails, is simply dealt with at the level of the UFO model [48] that MadGraph5 aMC@NLO has to import prior to generating a process (see appendix B.1 of ref. [27] for further details), and by including the relevant routines that perform the running of m b (µ). There is only one extra complication, due to the fact that the MS Yukawa introduces in the cross section an extra µ R dependence w.r.t. those taken explicitly into account in ref. [49], which are used by the code for the definition of scale-  and PDF-independent coefficients that are exploited for the a-posteriori computation of scale and PDF uncertainties by means of reweighting. Furthermore, such a dependence is different in the σ y 2 b and σ y b yt terms introduced in eq. (2.1), owing to the different powers of y b that appear in those terms. Although this complication will become recurrent in a mixed-coupling expansion scenario, at the moment it does not warrant a completely gen- eral and automated solution. Therefore, we have treated bbH production as a special case, by integrating separately the σ y 2 b and σ y b yt terms, which necessitate loop-content filtering operations (see sect. 2.4.2 of ref. [27]) in order to exclude, or to include only, top-quark loops in the virtuals. For each of these two terms, we have manually performed the modifications in the definition of the coefficients, mentioned above, that serve to compute the theoretical systematics. Apart from these manipulations, the generation and subsequent computation of the bbH 4FS cross section proceed exactly with the same general steps as those described in ref. [27], namely: MG5 aMC> import model loop sm MSbar yb MG5 aMC> generate p p > h b b~[QCD] followed by the standard output and launch commands, and where loop sm MSbar yb is the name of the UFO model that includes the appropriate UV counterterms for the renormalisation of the bottom Yukawa in the MS scheme.
Given that the few manual operations mentioned above are necessary on top of the commands just listed, the user interested in the simulation of bbH production with Mad-Graph5 aMC@NLO is strongly encouraged to contact us.

Five-flavour scheme
Consistently with the general discussion given in sect. 1, our 5FS results are obtained by setting m b = 0 (while keeping the bottom Yukawa finite). Sample Feynman diagrams that contribute to the cross section in this scheme are displayed in fig. 4a (at the LO, O(α 0 S )) and fig. 4b-4d (at the NLO, O(α S )); the cross section is proportional to y 2 b , and the y b y t term is absent at the α S order at which we are working. Analogously to what has been done in the 4FS, the bottom Yukawa is renormalised in the MS scheme (see sect. 2.1). The general comments made before that concern the generation of the process apply to the 5FS case as well. At variance with the 4FS calculation, however, the imported model must have the b quark mass set equal to zero (except in the Yukawa coupling). The relevant MadGraph5 aMC@NLO commands are thus: MG5 aMC> import model loop sm MSbar yb-no b mass MG5 aMC> define p = g u d s c b u~d~s~c~bM G5 aMC> generate p p > h [QCD] where the second line explicitly instructs the code to consider the b quark as part of the proton. Note that since the imported model is the Standard Model, and not the effective theory which features an effective ggH vertex, the generate command will indeed result in creating the 5FS bbH cross section we are interested into.

Phenomenological results
In this section we present several differential distributions that we reconstruct from finalstate particles in bbH production at the 13 TeV LHC. Although we work in the SM, our predictions are directly applicable to bbφ production (with a neutral φ = h, H, A) in a 2HDM-type extension of the SM, by an appropriate rescaling of the bottom Yukawa; in the case of the MSSM, this has been verified [50,51] to be an excellent approximation of the complete result.
Our reference predictions are (N)LO+PS-accurate; where appropriate, we also show f(N)LO and (N)NLO+(N)NLL results. For simulations matched to parton showers, we employ both Herwig++ [52,53] and Pythia8 [54]; for further information on the calculation of the MC counterterms relevant to the MC@NLO method, see refs. [23,55]. Higgs decays have not been considered; in particular, the contents of all jets in the events are solely due either to hard-process particles, or to radiation off those particles.

Input parameters
The central values of the renormalisation (µ R ) and factorisation (µ F ) scales for all Mad-Graph5 aMC@NLO runs ((N)LO+PS and f(N)LO) have been taken equal to the reference scale: where the sum runs over all final-state particles at the hard-process level; this is in keeping with the findings of refs. [19][20][21][22]. The theoretical uncertainties due to the µ R and µ F dependencies have been evaluated by varying these scales independently in the range: The calculation of this theory systematics does not entail any independent runs, and is performed by means of the reweighting technique introduced in ref. [49], with the bbHspecific upgrade discussed in sect. 2.1. In the case of the analytically-resummed cross sections, we have set the reference scale equal to: The so-called resummation scale (Q res ) in analytic transverse momentum resummation plays the role of a matching scale between the low-and high-p T regions. In order to optimise the high-p T matching of the resummed to the fixed-order cross section, this scale should be set equal to about half of the factorisation scale [26], and therefore we choose Q res = m H /8 as our default value; note that at the NNLO+NNLL, the choice of the resummation scale has hardly any impact in the small-p T region.
We have adopted the MSTW2008 PDF set [56], with its associated α S value, in its four-or five-flavour variant in agreement with the 4FS or 5FS computation being carried out. Jets are reconstructed with the anti-k T algorithm [59], as implemented in FastJet [60], with a jet radius of R = 0.5, and subject to the condition p T (j) ≥ 25 GeV. In the case of (N)LO+PS simulations, jets are made up of hadrons; b-jets (defined to be jets that contain at least one B hadron (at (N)LO+PS) or b quark (at f(N)LO)) are kept only if they fulfill the extra condition |η(j b )| ≤ 2.5. We have not generated underlying events.

Four-flavour scheme results
In this section we present 4FS results for total rates, possibly within cuts, and for differential distributions constructed with the four-momenta of the Higgs, b-jets, and B hadrons or b quarks. Before giving definite predictions for those quantities, however, there is a general issue that we would like to address, directly related to the novelty of being able to perform 4FS computations at the NLO+PS accuracy. This issue stems from a general characteristics of bbH production mentioned in sect. 1, namely the fact that the optimal values for the hard scales that enter the calculation appear to be significantly smaller than the hardness of the process would suggest, which is what has led us to the setting of eq. (3.1). When one considers parton showers, another hard scale becomes relevant, which loosely speaking can be identified with the largest hardness accessible to the shower; let us denote this scale by Q sh . It is the MC that determines, event-by-event, the value of Q sh , by choosing it so as to maximise the kinematic population of the phase-space due to shower radiation, without overstretching the approximations upon which the MC is based. The latter condition typically implies that where the average is taken over all generated events. In the context of the MC@NLO method, the condition of eq. (3.4) is actually not so crucial: essentially, MC hard radiation is subtracted, and replaced by that of matrix element origin, and in this way its impact on physical observables is suppressed. It is however not identically equal to zero, because the subtraction works within approximation, i.e. at the NLO; hence, MC hard radiation formally of NNLO and beyond can still contribute to the cross section. In order to assess this higher-order systematics of MC origin (which in MC@NLO is tantamount to the matching systematics), in MadGraph5 aMC@NLO one is given the possibility of setting the value 3 of Q sh ; this value is actually picked up at random in a user-defined range: so as to avoid possible numerical inaccuracies due to the presence of sharp thresholds; more details can be found in sect. 2.4.4 of ref. [27] (see in particular eq. (2.113) and the related discussion). In eq. (3.5), s 0 is the Born-level partonic c.m. energy squared, and α, f 1 , and f 2 are numerical constants whose defaults are 1, 0.1, and 1 respectively. The way in which Q sh is generated results in a distribution peaked at values slightly larger than α(f 1 + f 2 ) s 0 /2. The essence of the MC@NLO method is such that, in practice, virtually all processes studied so far exhibit a modest systematics due to the parameters that control Q sh . It is clear by construction that there is a rather direct relationship between Q sh in MC@NLO, and Q res in analytical resummation. Therefore, given the fact that bbH production in the 4FS is a chief example for the condition of eq. (3.4) not to be fulfilled, and that for this process one tends to use small values of Q res , it is interesting to investigate the sensitivity of NLO+PS predictions to the choices of the parameters that appear in eq. (3.5). Thanks to the redundancy of the latter, we shall limit ourselves here to studying the dependence on α, by setting α = 1, 1/2, 1/4. We shall consider an observable which, in NLO+PS computations, is maximally sensitive to the matrix-element vs parton-shower interplay, namely the transverse momentum of the Born-level "system" (p syst T ). The latter is unambiguosly defined only in fixed-order calculation, where its four momentum is the sum of the four-momenta of the Higgs, the b, and theb quark. In the case of NLO+PS simulations, we use the sum of the four-momenta of the Higgs and of the two hardest B hadrons; no final-state cuts are applied.
The results are presented in fig. 5, for both Herwig++ (left panel) and Pythia8 (right panel); the fNLO prediction is shown as well. The insets display the ratio of the NLO+PS results over the fNLO one. Only the y 2 b terms (σ y 2 b in eq. (2.1)) are considered here. We remind the reader that at sufficiently large transverse momentum the NLO+PS results (obtained with the MC@NLO method) will coincide by construction with the fNLO result (in the case of p syst T , up to small effects due to the fact that the hadron-and parton-level observables are not exactly the same). The main message one derives from fig. 5 is that shower (i.e. resummation) effects, measured by the distance between the NLO+PS and fNLO predictions, extend much farther than one would naively expect if "large" values of α are chosen. Furthermore, the dependence on α at large p syst T is extremely significant which, as explained above, is exceedingly rare for MC@NLO results. On the other hand, the NLO+PS curves do behave as expected: their shapes show no dependence on α at small p syst T , and their total integrals are equal at the level of the statistical integration error (i.e. 0.5%), and equal to the fNLO rate. These observations apply to both Herwig++ and Pythia8, which follow a rather similar pattern. What one sees, thus, is an extremely large matching systematics in certain corners of the phase space. On the one hand, this effect is enhanced by the fact that the tail of the p syst T distribution is extremely steep. On the other, the dynamics of the process is such that the condition of eq. (3.4) appears to be significantly violated. In particular, with α = 1 the distribution of Q sh peaks at around 180 GeV; this value decreases to 90 GeV and 45 GeV when setting α = 1/2 and α = 1/4 respectively. Therefore, with α = 1/4 one induces shower scales which are closer to the values taken by the hard scales, defined as in eq. (3.1).
At this point, there are two things which must be stressed. Firstly, by choosing α = 1 and µ as in eq. (3.1) one does not really introduce large logarithms in the computation. Rather, if the "small" value of µ is dictated by arguments of collinear dominance, the same arguments appear to suggest that, since the MCs are based on a collinear approximation, small values of Q sh have to be preferred. If these values are natural for bbH production, by using them one makes sure that the MCs radiate mostly away from the hard regions. An indirect confirmation of this can be seen in fig. 5: the smaller α, the closer Herwig++ and Pythia8 are to each other. Secondly, although from what has been said above a coherent picture emerges, one has to bear in mind that the entire discussion stems from the very large theoretical uncertainties that affect bbH production even at the NLO. So while this justifies, to a certain extent, the practice of finding optimal scale choices, it does not allow one to ignore the existence of large systematics, which might prove to be crucial in the comparisons with data.
Owing to the arguments above, α = 1/4 is our default choice in 4FS (N)LO+PS simulations, whose results we are now going to present. It is interesting, and reassuring for the self-consistency of the theoretical description of bbH production, that the hard scales relevant to our novel NLO+PS calculation appear to follow the same pattern as those relevant to other approaches to bbH production. In this paper, we shall refrain from quoting an uncertainty associated with the variation of α.
We start by reporting, in table 1, the results for the total cross sections, both fully inclusive and within cuts. As far as the latter are concerned, we have considered three possibilities: the requirement that there be at least one or two b-jet(s) (see sect. 3.1 for the jet-finding parameters), and that the transverse momentum of the Higgs be larger than 100 GeV (boosted scenario). In the case of the fully inclusive cross section, we have also computed it by setting µ R = µ F = m H /4, for ease of comparison with the results in the literature, with which we find agreement at the level of the numerical integration errors. For the cross sections within cuts, since they depend on the kinematics of final-state objects, we report the results obtained both with Pythia8 and with Herwig++. The fractional scale uncertainties that we show in table 1 are computed by varying the scales as indicated in eq. (3.2); since they are largely MC-independent, we give them only in the case of the Pythia8 simulations. We present separately the results for the y 2 b and y b y t terms. The conclusions that can be drawn from the table are the following.
1. The inclusion of NLO corrections reduces in a very significant manner the scale dependence of the y 2 b terms w.r.t. that at the LO; the residual scale dependence is however still large (∼ 20%).
2. The effect on the fully-inclusive y 2 b cross section of the NLO corrections is moderate (∼ 10%) and negative (i.e. the K factor is slightly smaller than one).
3. There are large differences (∼ 15%) between the cross sections computed with dynamic or fixed hard scales, of the same order as the NLO scale uncertainties.
4. The contribution of the y b y t term is negative and non negligible (∼ 10%), although within the scale uncertainties of the y 2 b NLO results. Its scale dependence is larger than that of the y 2 b term (as expected, owing to the absence of a y b y t contribution at the LO).

In the cases of the cross sections within cuts, the inclusion of NLO corrections im-
proves the agreement between the predictions of the two MCs w.r.t. that at the LO.
6. The effects of the cuts are significant: the cross section is reduced by a factor larger than three when requiring at least one b-jet, and by a further factor of about ten when a second b-jet is tagged. The boosted-Higgs case is similar to the two-jet one.
We now turn to discussing differential observables. Unless stated otherwise, we shall limit ourselves to presenting the results relevant to the y 2 b terms since, apart from these being dominant, it turns out that the y b y t contribution is fairly flat for most observables; there is one striking exception, which we shall discuss in details later. The figures of this section and of sect. 3.3 are generally organised according to the following pattern. There is a main frame, where the relevant predictions (e.g. NLO+PS, fNLO, and so forth) are shown with their absolute normalisation, and as cross section per bin (namely, the sum of the contents of the bins is equal to the total cross section, possibly within cuts). In an inset we display the bin-by-bin ratio of all the histograms which appear in the main frame over one of them, chosen as a reference. Finally, in a second inset the bands that represent the fractional scale dependence are given: they are computed by taking the bin-by-bin ratios of the maximum and the minimum (obtained according to eq. (3.2)) of a given simulation over the same central prediction that has been used as a reference for the ratios of the first inset. In this paper we have ignored the PDF systematics; however, we stress that it can be included at no extra CPU cost, using the same reweighting procedure [49] employed here for the scale uncertainty.
We begin in fig. 6 with the transverse momentum of the Higgs. In the left panel we present the NLO+PS Herwig++ result (black solid; this is the reference histogram), its LO+PS counterpart (red dotted), and the fNLO (blue dot-dashed) and fLO (green dash-double-dotted) predictions. In the upper inset, on top of the ratios of the curves that appear in the main frame we also display the ratio of the NLO+PS Pythia8 result (orange solid with full boxes) over the Herwig++ one. In the region p T (H) 20 GeV, all results are within 25% of each other. The agreement between the NLO+PS Herwig++ and Pythia8 predictions is excellent (they are essentially identical up to statistics); the fNLO result is slightly harder than its showered counterparts, and it gets closer to them with increasing p T (H). The LO+PS Herwig++ prediction is within 10% of the NLO+PS one in the whole p T (H) range. Pythia8 differs visibly from Herwig++ at the NLO+PS only in the small p T (H) region, where the difference is however 20% at most; on the other hand, in that region it is the fixed-order results which display the largest discrepancies w.r.t. our Herwig++ NLO+PS reference curve. From the lower inset one sees the very significant reduction of the scale dependence at the NLO w.r.t. the one at the LO, analogous to that already observed in table 1 for total cross sections. The uncertainty bands of the NLO predictions almost completely overlap, and they are by far and large contained within those relevant at the LO.
The Herwig++ and Pythia8 NLO+PS results, and the fNLO one, are also shown in the right panel of fig. 6, with the same patterns as those employed in the left panel of that figure; on top of these, we also present the Herwig++ (red dotted) and Pythia8 (green dash-double-dotted with open boxes) NLO+PS predictions obtained with α = 1. Although the effects due to the change of α are not as large as those affecting p syst T (see fig. 5), they are very significant in the large-p T (H) region, where the two NLO+PS results obtained with the same MC have non-overlapping scale-uncertainty bands (the more so for Pythia8). Furthermore, while when α = 1/4 Herwig++ and Pythia8 agree perfectly with each other, this is not the case when α = 1 (although their uncertainty bands do overlap). Note that, as expected, the value of α does not affect the shape of the smallp T (H) distribution for a given MC. Overall, these results confirm the findings discussed above for p syst  are presented as black solid (NLO+PS) and red dotted (LO+PS) histograms, while those for b-jets are the blue dot-dashed (NLO+PS) and green dash-double-dotted (LO+PS) histograms. In the insets, the red dotted (green dash-double-dotted) histograms display the ratios of the LO+PS over the NLO+PS predictions for B hadrons (b-jets). In the case of Herwig++, the inclusion of the NLO corrections is hardly visible, except at the threshold of the B-hadron distributions, where the differences account for the smaller-than-one K factor reported in table 1 for total cross sections (the NLO+PS results being smaller than the LO+PS ones). The b-jet distributions are essentially identical, in shapes and rates, at the NLO+PS as at the LO+PS. In the case of Pythia8, one observes similar effects as for Herwig++ at the thresholds of the B-hadron distributions, but also a softening of the high-p T spectra at the NLO+PS w.r.t. the LO+PS, softening which is more pronounced in the case of the second-hardest B hadron. As far as b-jets are concerned, the NLO+PS and LO+PS results are quite similar in shape, but the rates of the former are about 10% and 20% lower than those of the latter for the hardest and second-hardest jet, respectively. Still, for each observable (a fortiori for Herwig++) the NLO+PS uncertainty band lies completely within that of the corresponding LO+PS prediction. The comparison of the two MCs for a given observable is best read from the upper insets, where we report the ratios Pythia8/Herwig++ ( fig. 7) and Herwig++/Pythia8 ( fig. 8) as solid orange (overlayed with boxes) and solid magenta (overlayed with stars) histograms for B hadrons and b-jets, respectively; we limit ourselves to presenting these ratios at the NLO+PS accuracy. For p T 's larger than about 20 GeV Herwig++ and Pythia8 agree quite well with each other in terms of shapes; rate-wise, Pythia8's are larger than Herwig++'s, with the largest differences (about 30%) occurring for the second-hardest B-hadron, and the smallest (about 10%) for b-jets. Significant differences (i.e. outside of the uncertainty bands) are seen only in the p T ∼ 0 region for B hadrons, where Pythia8 can be up to a factor of two smaller than Herwig++ in the first bin of the hardest-B hadron distribution.
The vast majority of the differential observables which we have investigated follow the pattern that underpins figs. 6-8: the various approximations are consistent within their associated theoretical uncertainties, and there are small corrections due to NLO and/or parton-shower effects and an overall agreement between the Herwig++ and Pythia8 predictions. A couple of exceptions are however notable, which we are now going to address explicitly.
In fig. 9 we present various predictions for the invariant mass of the two hardest B hadrons at the NLO+PS (black solid) and LO+PS (red dotted) accuracy, and that of the two b quarks at the fNLO (blue dot-dashed) and fLO (green dash-double-dotted) accuracy; no cuts are applied, and both quantities have been labelled m BB for simplicity. The left and right panels display the Herwig++ and Pythia8 results respectively. In the upper insets, on top of showing the ratios of the histograms that appear in the main frames over the relevant reference curves (Herwig++ NLO+PS and Pythia8 NLO+PS in the left and right panels respectively), we also show the ratios Pythia8/Herwig++ (left panel) and Herwig++/Pythia8 (right panel) at the NLO+PS (solid orange overlayed with boxes) and at the LO+PS (solid magenta overlayed with stars) accuracy. The most obvious feature of fig. 9 is the behaviour at threshold, where the differences between the two MCs are extremely large, and to a good extent independent of the perturbative order considered. The LO+PS Herwig++ prediction is quite far from the fLO result, which in turn is in reasonable agreement with Pythia8's. One should be careful not to interpret this fact as an evidence for the Pythia8 curve to be more realistic: the threshold region is where one expects to be more sensitive to the details of the b → B hadronisation, and thus where the differences between the invariant mass of the B hadrons and that of the b quarks are largest. The inclusion of NLO effects, regardless of the presence of parton showers, hardens the spectra in a very significant manner; owing to the shape of the LO+PS Herwig++ result, the shift in the peak position when going from LO to NLO is particularly large in the case of that MC. A remarkable fact is the similarity of the pattern LO+PS→NLO+PS in Pythia8 and Herwig++, in spite of the differences of the individual results at a given order: this can be seen clearly from the red dotted histograms in the upper insets of the two panels. Away from the threshold region, the two MCs are quite close to each other (and close to the fixed-order result, too) at the NLO; the agreement is worse, but still amply within theoretical uncertainties, at the LO, where however the fLO prediction is clearly farther away from the Herwig++'s than from the Pythia8's. Overall, Pythia8 results are harder than those of Herwig++, but the difference in hardness descreases with the perturbative order. We finally point out that the differences in the threshold region we have discussed above are greatly reduced or absent when one requires the presence of at least one large-p T b quark (e.g. in the invariant mass of two b-jets, or when cutting away the low-p T (B) region).
The second interesting, if less spectacular, case of significant differences among predictions is presented in fig. 10, where we show the distance in the η−ϕ plane between the two hardest B hadrons (at (N)LO+PS) or the two b quarks (at f(N)LO); both are denoted Figure 10: Same as in fig. 9, for the distance in the η−ϕ plane. by ∆R BB for simplicity. As in the case of m BB previously discussed, no cuts are applied. The layout of fig. 10 is the same as that of fig. 9. The salient feature of the present results is to be found at large ∆R BB (which corresponds to a large separation in pseudorapidity between the two B hadrons or b quarks). In that region, both MCs tend to be lower than the corresponding fixed-order result; however, in the case of Pythia8 the decrease of the cross section is dramatic (even at the NLO, where it is nonetheless smaller than at the LO), while it is modest (and within scale uncertainty) for Herwig++. Over the whole spectrum, NLO corrections have a large impact, depleting the cross section at small ∆R BB and enhancing it at large ∆R BB . This pattern is independent of whether parton showers are included; note that in some regions the LO and NLO uncertainty bands do not overlap. Up to ∆R BB ∼ 6 and ∆R BB ∼ 7.5 (at the LO and NLO, respectively), the mutual agreement between MCs and with fixed-order results is quite good, but it rapidly degrades above those values as was already discussed; for all ∆R BB , however, the agreement at the NLO is better than that at the LO. We finally remark that, in terms of the comparison of the MC predictions with the fixed-order ones, the situation at large ∆R BB is the opposite as that for m BB at threshold, with Pythia8 (N)LO+PS being farther away from f(N)LO than Herwig++ for the former observable, and closer for the latter. While, as was stressed before, the level of agreement with f(N)LO results cannot be used as a discriminant for MCs, this fact underlines the different characteristics of different MCs, and the necessity of a conservative estimate of MC systematics.
We now turn to discussing the impact of the y b y t term (i.e. the σ y b yt contribution to the NLO cross section in eq. (2.1)). We have investigated a rather large number of observables, and generally found such a term to be flat at the level of 5%-20% (this fraction is measured by evaluating the ratio of the y b y t over the y 2 b contribution, and by taking its distance from a horizontal line); we remind the reader (see table 1) that the y b y t term gives a negative contribution of the order of 10% at the level of total rates. Therefore, although there is no reason to neglect the y b y t contribution, an overall rescaling is a decent approximation (also in view of the large scale uncertainties that affect the y 2 b bbH cross section) for most observables. There are two counterexamples, one of which is particularly spectacular, which we present below; perhaps not surprisingly, they coincide with the m BB and ∆R BB observables which we have just discussed.
The results for m BB are displayed in fig. 11, at the NLO+PS accuracy with Herwig++. The black solid histogram is the result for σ y 2 b , and the red dotted histogram the result for σ y 2 b + σ y b yt . The left panel presents the inclusive case, while in the right panel we have required the presence of at least one b-jet. We start by commenting the inclusive case: as one can see, the two predictions are in very good agreement at large m BB (up to a rescaling), and in violent disagreement close to threshold, where the striking feature is a sharp peak which originates from the y b y t contribution. Given the situation of fig. 9, and in particular the large differences between Herwig++ and Pythia8 there, in the upper inset we report the ratio (σ y 2 b + σ y b yt )/σ y 2 b not only for Herwig++ (i.e. the ratio of the two curves displayed in the main frame; red dotted), but also for Pythia8 (blue dot-dashed). As one can see, the two ratios are in rather good agreement with each other; in other words, the pattern of the peak vs no-peak structure close to the threshold of m BB is essentially MC-independent, and is thus purely of matrix-element origin. It is in fact straightforward Figure 12: Same as in fig. 11, for the distance in the η−ϕ plane.
to connect this behaviour to the topology of the one-loop graphs that contribute to the y b y t term (see for example fig. 2a, 2b, and 2e): the vast majority of them feature a g → bb splitting, which is naturally enhanced at small m BB (in the case of the y 2 b contribution, only a small fraction of diagrams contain such a splitting). The enhancement due to the g → bb splitting is easily countered, for example by requiring the presence in the event of a b-jet: the result is shown in the right panel of fig. 11, where the sharp peak at small m BB does not appear any longer; the effect of the y b y t terms is now much milder (and the relative behaviour of Herwig++ and Pythia8 is again almost identical). We have verified that, in a completely analogous manner, the y b y t contribution is quite flat also in the case of the invariant mass of the two hardest b-jets.
The case of ∆R BB is presented in fig. 12, which has the same layout as fig. 11. The impact of the y b y t contribution is not as outstanding for the present observable as for m BB , but it is still clearly visible, in that it causes the σ y 2 b + σ y b yt cross section to be larger than the σ y 2 b one in the small ∆R BB region, in the inclusive case (see the left panel of fig. 12). As already for m BB , the requirement that there be at least one b-jet changes the picture (see the right panel of fig. 12); still, the y b y t contribution is less flat than for the vast majority of the other observables, and is at the border of the σ y 2 b scale-dependence band.

Four-and five-flavor scheme comparison
In this section we present selected 5FS predictions, and compare them directly to their 4FS counterparts, many of which have been already shown in sect. 3.2. We shall mostly work at the NLO+PS level, although occasionally we shall use LO+PS results as well; furthermore, where appropriate we shall also compare our predictions to those of (N)NLO+(N)NLL analytical resummations [26]. Since the y b y t term vanishes in the 5FS up to O(α 2 S ), we shall only consider the y 2 b contribution throughout this section. We also point out that, in the case of a 5FS computation, one has s 0 = m 2 H (see eq. (3.5)). Therefore, at variance with the case of the 4FS, Q sh naturally assumes values rather smaller than the Higgs mass given the default f i parameters. Thus, in the present case the default choice α = 1 is perfectly adequate; we have in any case verified that the dependence on α is moderate (O(10%), and affecting mostly the intermediate-p T region), and much smaller than for 4FS results.
We start by considering the Higgs transverse momentum. In the left panel of fig. 13 the NLO+PS Herwig++ and Pythia8 5FS predictions (blue dot-dashed and green dashdouble-dotted respectively) are compared with the results of the analytical resummation at the NLO+NLL (red dotted) and NNLO+NNLL (black solid) accuracy; the latter result is the reference curve, used as the denominator in the computation of the ratios which appear in the insets. The analytically resummed results, in this plot and in all those which will follow, have been rescaled by the bin size for a direct comparison of both rate and shape with the NLO+PS predictions (since the latter are always in the form of cross sections per bin, while the former are returned by the code of ref. [26] as differential distributions); their associated uncertainty bands (which in the case of the NNLO+NNLL result are smaller than those relevant to any other simulation considered here) include an independent variation of Q res . In the small-p T region, there is a good agreement among the two NLO+PS and the NLO+NLL results, while the NNLO+NNLL one is visibly lower. The peaks of all curves lie within 5 GeV of each other, that of Herwig++ (NNLO+NNLL) being the lowest (highest). Starting from about p T (H) ∼ 60 GeV, the NLO+PS results are closer to the NNLO+NNLL curve than to the NLO+NLL one (which however is within the scale uncertainty bands of the former), the agreement being particularly good in the case of Pythia8. Note that at large p T the NNLO and NLO predictions are quite close to each other; this is analogous to what has been observed in ref. [42] for the transverse momentum of the hardest jet, and is a consequence of using µ R = µ F = m T (H)/4. In the right panel of fig. 13 we compare the most accurate 5FS prediction, namely the analyticallyresummed NNLO+NNLL, with the NLO+PS Herwig++ and Pythia8 ones in the fourflavour scheme, which have already appeared in fig. 6. As we know from that figure, the agreement between the two NLO+PS results is excellent; what one sees from fig. 13 is that these NLO+PS predictions also agree rather well with the NNLO+NNLL one (and in an excellent manner shape-wise), except when very close to p T (H) = 0, with fully overlapping uncertainty bands. In fact, up to an overall rescaling of the 4FS NLO+PS curves, one may say that the agreement with the NNLO+NNLL result is better for the 4FS NLO+PS results than for the 5FS NLO+PS ones. However, one must bear in mind that the NLO+PS Figure 15: Rapidity of the hardest (left panel) and second-hardest (right panel) B hadron, in the 4FS and 5FS at the NLO+PS accuracy, as predicted by Herwig++ and Pythia8. All histograms have been normalised so that their integrals are equal to one. matching systematics (see eq. (3.5)) is much larger in the 4FS than in the 5FS. Conversely, note that widths of the NLO+PS 5FS uncertainty bands are larger than those relevant to the 4FS for p T (H) 80 GeV, because from the perturbative viewpoint that kinematic region is effectively described at the LO in the 5FS.
The p T (H) distribution is severely affected by the requirement that there be at least one b-jet in the final state. In fig. 14 we present the relevant results, obtained at the NLO+PS accuracy in the 4FS (black solid (Herwig++) and red dotted (Pythia8)) and in the 5FS (blue dot-dashed (Herwig++) and green dash-double-dotted (Pythia8)). For p T (H) 50 GeV all predictions are within 30% of each other, the agreement among the two 4FS results and the 5FS Herwig++ one being particularly good. Below 50 GeV more significant deviations (especially in shape) start to appear, which increase with decreasing p T (H). The pattern of the comparison between the two MCs is the same in the two schemes: namely, the Pythia8 cross section is larger than the corresponding Herwig++ one. The two 4FS results are larger than the two 5FS ones. Below p T (H) ∼ 10 GeV, the uncertainty bands of the 4FS results do not overlap any longer with those of the 5FS ones; within a given scheme, the bands do overlap, but the central predictions show differences at the level of 30% and 20% in the 4FS and 5FS respectively. In conclusion, although the overall agreement between the 4FS and 5FS results is reasonable, shape-wise visible discrepancies do appear, which would thus be interesting to investigate using data, especially in view of the fact that theoretically, for an observable that features a tagged b-jet, the 4FS is expected to be superior to the 5FS. In the context of MC simulations the presence of a massless b quark in the initial state at the matrix-element level poses a non-trivial problem; apart from the necessity of evolving it backwards in a way that matches the flavour content of the incoming hadron, one always faces a kinematic constraint, imposed by the fact that eventually the b quark will have to appear in the final state, with a mass of around 5 GeV. This problem, which is essentially process-independent, is particularly severe in the case of HERWIG6 [61,62], where it leads to strangely-looking distributions, with longitudinal observables being particularly affected. An explicit example is given in fig. 14 of ref. [23], where 5FS single-top production has been considered; from that figure, one can also see how this issue has been addressed by Herwig++, which (almost completely) rectifies the behaviour of HERWIG6. It is important to stress that it is the region p T ∼ 0 which is relevant here: as soon as one imposes a realistic tagging condition (by requiring a minimum p T in order to observe a B hadron), or is sufficiently inclusive, the problem above becomes essentially irrelevant. Still, from the theoretical point of view this small-p T behaviour of b-flavoured hadrons is less than satisfactory.
With this in mind, and given the improvement of Herwig++ w.r.t. HERWIG6 in the case of single-top production, we have considered the rapidity distributions of the hardest and second-hardest B hadrons. The (normalised) results are displayed in the left and right panel, respectively, of fig. 15, where we show both the 4FS predictions (black solid (Herwig++) and red dotted (Pythia8)) and the 5FS ones (blue dot-dashed (Herwig++) and green dash-double-dotted (Pythia8)). The prominent feature of these plots is the behaviour of the Herwig++ 5FS results, which are vastly broader than all of the other Figure 17: Comparison between the 4FS and 5FS results for the invariant mass of (left panel) and the η−ϕ distance between (right panel) the two hardest B hadrons. three curves. Therefore, although no features appear such as those mentioned above for HERWIG6, it is likely that Herwig++ still tends to produce B hadrons too close to the beam line when simulations are performed in the 5FS. As far as the other histograms are concerned, the agreement between the two Pythia8 predictions is extremely good for rapidities as large as about 3. Beyond this value, the 5FS result is broader than the 4FS one, and is actually rather close to the Herwig++ 4FS result in the case of the hardest B hadron (not so for the second-hardest B hadron, for which at large rapidities the two 4FS predictions agree at the level of 20%).
In order to further the study of the behaviour of Herwig++ for these rapidities, in fig. 16 we present the comparison between the NLO+PS and LO+PS predictions, in both the 4FS and 5FS. As one can see, the effect of the NLO corrections is fairly modest: the 5FS NLO results are slightly more central than their LO counterparts, and the opposite happens in the 4FS. This implies that the two NLO+PS predictions are marginally closer to each other than the two LO+PS ones, but this fact is quite irrelevant given the vastly different shapes one obtains in the two schemes.
The observations above have bearings on the predictions for the two quantities which have been already addressed at length in sect. 3.2, namely the invariant mass of (m BB ), and the η − ϕ distance between (∆R BB ), the two hardest B hadrons. We display these observables in fig. 17, by comparing the 4FS predictions already presented in fig. 9 and fig. 10 to their 5FS counterparts. The most evident feature in the two panels of fig. 17 is the very large Herwig++ 5FS cross section at large m BB and ∆R BB , which suggests again that Herwig++ in this scheme tends to produce B hadrons fairly close to the beam line, and in opposite directions. The other three results are much closer to each other, although large differences remain; one may note how the 4FS Herwig++ predictions have shapes rather similar to the 5FS Pythia8 ones. There is no fundamental reason why this should be so; in particular, one should bear in mind that for these observables one expects the 4FS to be significantly more reliable than the 5FS. Still, our results indicate that the underlying MCs, even in the context of matched simulations at the NLO, play a non-negligible role. A thorough comparison with data, for this and other b-initiated processes, will certainly be beneficial for a better understanding of these issues, and for improving the tuning of the relevant long-distance parameters.

Conclusions
In this paper we have studied the hadroproduction of a Higgs boson in association with b quarks, and presented for the first time NLO QCD predictions matched to parton showers. We have worked in the MadGraph5 aMC@NLO framework, endowed with the capability, which was unavailable in the code prior to this study, of treating the renormalisation of the bottom Yukawa coupling in a fully flexible manner, and in particular in the MS scheme. We have provided results in both the four-flavour and the five-flavour schemes, which we have compared to each other and, in the case of the Higgs transverse momentum, with analytically-resummed (N)NLO+(N)NLL results as well. We have also considered the O(y b y t α 3 S ) term in the 4FS, which might be viewed as the leading contribution to the interference between the bbH and gluon-fusion channels.
The key feature of the predictions we report is their being fully differential, independently of whether they are matched with parton showers. We have documented this by discussing the cases of several observables that are exclusive in the Higgs and/or in the b-flavoured hadron or b-quark momenta, and which are notable for their characteristics. Although such observables represent only a limited sample of what can be obtained by means of MadGraph5 aMC@NLO, they extend in a very significant manner the knowledge of differential quantities which was available in the literature so far. One aspect of our results that is particularly worth stressing is that their associated scale (and PDF, which for simplicity we have refrained from investigating here) uncertainties can be computed at the same time as the central results without any noticeable CPU cost. This is particularly important in the case of bbH production, in view of the large theoretical systematics which affect this process, and which must thus be carefully taken into account, as we have done for all predictions given in sect. 3.
Because of the fully-exclusive nature of our computations, we regard the 4FS results, especially thanks to the matching to parton showers, as generally superior to the 5FS ones, and we believe that they should constitute the default choice for any realistic physics simulation. Having said that, for observables that are fully inclusive in the degrees of freedom of the associated b quarks, or for which m b is negligible (e.g. at large p T 's) the differences between the four-and five-flavour schemes must be carefully assessed. From the phenomenology viewpoint, the main conclusions of this work are the following.
• We have found evidence that relatively small values for the shower scales in 4FS computations have to be preferred. This is in keeping with the by-now standard choices for the other hard scales that enter the process, and with analogous findings in the context of other types of calculations (i.e. fixed order in both the four-and five-flavour schemes, and analytical resummations).
• The impact of NLO QCD corrections is very significant, both in terms of reducing the scale uncertainties w.r.t. those that affect the LO results, and in the changes they induce in the shapes of some differential observables (in particular, the m BB and ∆R BB correlations). Even at the NLO, however, the perturbative systematics that affect bbH production are sizable.
• The O(y b y t α 3 S ) interference contribution reduces the inclusive rate by about 10%. It is generally flat in the phase space, except close to the threshold of the invariant mass of the hardest-B-hadron pair, where it gives a very prominent peak structure. Such a peak tends to disappear when increasing the minimum transverse momentum of the B hadrons, or when b jets are employed.
• The matching of NLO results with parton showers plays a very important role. On the one hand, there are observables for which fixed-order predictions are significantly different w.r.t. those after showers. On the other hand, there are cases where the Pythia8 and Herwig++ results show large discrepancies, owing to fundamental differences in the implementation of core physics characteristics (such as shower and hadronisation mechanisms); this is especially true for 5FS simulations. The MC systematics must thus be considered very carefully in an observable-by-observable manner.
• The NLO+PS predictions for the Higgs transverse momentum, inclusive in the degrees of freedom of the b quarks, compare remarkably well with the analyticallyresummed ones. In terms of shape, this is particularly true for the 4FS and the NNLO+NNLL results.