RIP $H b \bar b$: How other Higgs production modes conspire to kill a rare signal at the LHC

The hadroproduction of a Higgs boson in association with a bottom-quark pair ($H b \bar b$) is commonly considered as the key process for directly probing the Yukawa interaction between the Higgs boson and the bottom quark ($y_b$). However, in the Standard-Model (SM) this process is also known to suffer from very large irreducible backgrounds from other Higgs production channels, notably gluon-fusion ($gg$F). In this paper we calculate for the first time the so-called QCD and electroweak complete-NLO predictions for $H b \bar b$ production, using the four-flavour scheme. Our calculation shows that not only the $gg$F but also the $ZH$ and even the vector-boson-fusion channels are sizeable irreducible backgrounds. Moreover, we demonstrate that, at the LHC, the rate of these backgrounds is very large with respect to $H b \bar b$ production and in particular no suppression occurs at the differential level. Therefore, they survive typical analysis cuts. This fact further jeopardises the experimental measurement of SM $H b \bar b$ production at the LHC. Especially, unless the $y_b$ is significantly enlarged by new physics, even for beyond-the SM scenarios the direct determination of $y_b$ via this process seems to be hopeless at the LHC.


Introduction
Since the beginning of its operation the Large Hadron Collider (LHC) has disclosed an unprecedented amount of information on the nature of the Higgs boson. After having discovered this particle [1,2], the ATLAS and CMS collaborations have also observed all the four main production modes [1][2][3][4][5][6][7][8]: gluon-fusion (ggF), vector-boson-fusion (VBF), vector-boson associated production (V H), and top-quark pair associated production (ttH). Moreover, many of the decay channels, γγ, V V * , τ + τ − and bb have been observed, too. All these measurements has lead to a very clear statement: the properties of this particle are compatible with those of the Higgs boson of the Standard Model (SM) of elementary particles. Notably, the couplings of the Higgs boson have already been found to be compatible with the SM predictions at 15-20% accuracy for the case of the third-generation fermions and even at ∼ 5% accuracy for the W and Z bosons [9].
On the other hand, there is still large room for possible beyond-the-SM (BSM) effects. For instance, with the exception of the muon [10,11], couplings to the first two fermion generations are largely unconstrained. A similar situation is also present at the moment for the Higgs self coupling, for which only bounds of the order of 5-10 times the SM value have been obtained [12,13]. In addition, the extracted value for a specific Higgs coupling depends on the assumptions on (the structure of) other interactions, both in the Higgs sector as well as in the top-quark or electroweak (EW) sector. For instance, Higgs production processes at the LHC are always measured in conjunction with specific decay channels, including those still unconstrained. Therefore the extraction of the couplings in a given production mode is affected by those appearing in the decay channel, and vice versa. Moreover, branching ratios depend on the total decay width and in turn on all the other decay channels. For this reason, having independent measurements of various production processes can certainly improve the precision in the determination of a specific coupling, and meanwhile it allows to relax the underlying theoretical assumptions.
At the LHC, in the case of the bottom-quark Yukawa coupling (y b ), a direct sensitivity can be in principle achieved both via the H→bb decay or via the associate production of a Higgs boson together with a bb pair. While the H→bb decay has already been observed in conjunction with V H production [5,6,14,15], no dedicated SM analysis has ever been performed by the ATLAS and CMS collaboration in order to measure the Hbb production process. Indeed, this process has an inclusive cross section that is comparable to the one of ttH production (e.g. ∼ 0.5 pb at 13 TeV). However, at variance with ttH production, at least one b-jet has to be tagged in order to make it distinguishable from the inclusive Higgs-boson production process, whose production rate is ∼ 100 times larger than Hbb production alone. The problem is that tagging a b-jet dramatically reduces the cross section, but as said this is an unavoidable procedure for obtaining a possible direct sensitivity on the bottom-quark Yukawa coupling. Without tagging any b-jet, even bottom-quark loops in ggF, which induce an indirect sensitivity on y b , have a larger contribution: about −6% for the inclusive cross section and up to −10% for the Higgs boson at small transverse momentum [16][17][18][19][20][21][22][23].
In the past, Hbb production has nevertheless received a lot of attention from the theoretical community. Indeed, in BSM theories with an extended Higgs sector, such as the two-Higgsdoublet-model (2HDM) or the minimal-supersymmetric-SM (MSSM), the coupling of the Higgs boson to bottom quarks can be significantly enhanced. Moreover, in the context of higherorder QCD corrections, this process is particularly interesting also from a formal point of view. Featuring bottom quarks in the hard process, two different schemes can be adopted when performing perturbative calculations: the so-called four-flavour scheme (4FS) and five-flavour scheme (5FS). In the former, the bottom quark is considered as massive, while in the latter the bottom mass is set equal to zero.
In the 4FS, next-to-leading order (NLO) QCD predictions were firstly obtained in Refs. [24,25], then for MSSM-type couplings in Ref. [26], and including supersymmetry (SUSY) QCD corrections in Ref. [27,28]. In the 5FS, many more results are present in the literature, since the higher-order perturbative calculations are technically much easier, owing to the smaller number of final-state particles. Indeed, corrections up to next-to-next-to-next-to-leading order (N 3 LO) QCD accuracy [29][30][31][32] are available. Distributions at the parton level were obtained for H+b and H+jet production at NLO in Refs. [33,34]. Next-to-next-to-leading order (NNLO) accuracy has been then reached for jet rates in Ref. [35] and for fully differential distributions in Ref. [36]. The spectrum of the Higgs boson transverse-momentum was studied analytically up to O(α 2 s ) in Ref. [37], while resumming at NLO+NLL and NNLO+NNLL accuracy in Refs. [38] and [39], respectively. For what concerns NLO QCD predictions matched to parton shower effects (NLO+PS), both the 4FS and the 5FS cases were studied for the first time within MadGraph5_aMC@NLO [40] and subsequently also via the Powheg [41] and Sherpa [42] approaches. Finally, at the level of the total cross section, differences between results obtained in the two schemes have been studied in Refs. [43,44] and then the state-of-the-art 4FS and 5FS predictions have been combined in Refs. [45][46][47][48] and very recently in Ref. [49], which matches the 5FS prediction at N 3 LO QCD accuracy and the 4FS prediction at NLO QCD.
The difficulties in the extraction of y b via the measurement of Hbb production do not originate only from its very small cross section. In the 4FS, NLO QCD corrections to Hbb interfere with gluon-fusion Higgs production (at LO) with an extra emission of a bb pair (ggF+bb). This interference leads to a term proportional to y b y t , where y t is the top-quark Yukawa coupling [40]. Moreover, the LO diagrams of ggF+bb are infrared (IR) finite and lead to terms proportional to y 2 t , when squared. The y t y b term is non-negligible w.r.t. the term proportional to y 2 b originating from the "genuine" Hbb production, and especially the y 2 t term is much larger than the y 2 b contribution. Both the y t y b and y 2 t terms have been calculated at NLO QCD accuracy in Ref. [50]; if at least one b-jet is required, the y t y b term is ∼ −20% of the y 2 b term, while the y 2 t term is ∼ 5 times the y 2 b term. In this paper, we calculate the cross section of Higgs boson production in association with a bb pair, Hbb, at "complete-NLO" accuracy in the 4FS. In other words, we compute all the SM contributions from QCD and EW origin at the tree and one-loop level. The choice of the 4FS is driven by the fact that when EW corrections are taken into account, the SM relation between y b and m b cannot be ignored and therefore in the 5FS one must enforce y b = 0. Complete-NLO predictions take into account not only the LO (order α 2 s α), the NLO QCD (order α 3 s α) and NLO EW (order α 2 s α 2 ) corrections, but also all the possible terms of order α m s α n+1 with m, n ≥ 0 and m + n = 2, 3. Part of the NLO EW corrections (only the gluon-gluon initial-state contribution) has already been calculated also in Ref. [51]. By considering complete-NLO predictions, new topologies open up on top of the aforementioned ggF+bb one. Terms with n ≥ 2 include ZH production, with subsequent Z decay into a bb pair. Moreover, at order α s α 3 even VBF configurations with Z bosons arise. In this paper we show, for the first time, that the numerical impact of these additional contributions is sizeable and sometimes even dominate, though they are formally suppressed by the small EW coupling constant. Especially, we show that the suppression of their relative contributions via ad hoc cuts inevitably strongly reduces also the total rates.
In our study, we demonstrate that the idea of directly extracting y b from the measurement of Hbb at the LHC is substantially hopeless. The rates for this process are small and contaminated by terms that depend on y t and the HZZ coupling. Reducing this contamination implies also a strong reduction of the cross section of the term depending only on y b .
The aforementioned computation has been performed via the latest version of MadGraph5_-aMC@NLO [52], which is public and has been extended in order to be able to calculate NLO EW corrections, and in general complete-NLO predictions, also in the 4FS. 1 Since the results presented in this paper represent the first complete-NLO (and also NLO EW) computation performed in such a scheme, we will also discuss in the text the relevant technical aspects connected to the usage of the 4FS in NLO EW corrections.
The paper is organised as follows. In Sec. 2 we describe our calculation set-up. The technical improvements performed for calculating complete-NLO predictions in the 4FS via MadGraph5_aMC@NLO are documented in Sec. 2.1. In Sec. 2.2 we describe the different topologies entering our calculation at different perturbative orders and in Sec. 2.3 we discuss the problems related to the MS renormalisation of y b , when EW corrections are present. Numerical results are presented in Sec. 3. Input parameters are given in Sec. 3.1, while numerical results at the inclusive and differential level are presented and commented in details in Secs. 3.2 and 3.3, respectively. The main phenomenological result of our work, i.e. the fact that the idea of directly extracting y b from the measurement of Hbb at the LHC is substantially hopeless, is motivated in detail in Sec. 3.2 and further corroborated by the analysis at the differential level in Sec. 3.3. We give our conclusions and outlook in Sec. 4.
2 Complete-NLO predictions for Hbb production As mentioned in Sec. 1, in this work we present the complete-NLO predictions for Hbb production at the LHC, namely, we exactly take into account all the one-loop and real-emission corrections of QCD and EW origin. Expanding in powers of α s and α, the first non-vanishing contribution to the cross section of Hbb production is of O(α 2 s α), and it is induced by tree-level gg → Hbb and qq → Hbb diagrams. Complete-NLO predictions for Hbb production include all the O(α m s α n+1 ) contributions with m, n ≥ 0 and m + n = 2, 3. Following the notation already used in Refs. [52][53][54][55][56][57][58][59][60], the different contributions to any differential or inclusive cross section Σ can be denoted as: The "standard" LO contribution is in our notation the LO 1 , while the term "LO" is rather used for denoting the sum of all LO i terms, consistently with Eq. (1). Similarly, according to Eq. (2), the term NLO is used to denote the sum of all NLO i terms. Thus, the quantity LO + NLO corresponds to the complete-NLO predictions. Further definitions will be given also later in the text in Eqs. (5)-(9); they constitute the quantities entering the phenomenological discussion of Sec. 3.
The calculation has been performed via the latest version of the public code MadGraph5-_aMC@NLO [52] (MG5_aMC hence-force), which enables the user to compute predictions at NLO EW and complete-NLO accuracy also in the 4FS, for a generic process in the SM or in any model for which the necessary counterterms are known. No ad hoc customisation of the code has been put in place in order to perform this specific calculation for the Hbb final state.
Technical details about the evaluation of loops and the ultraviolet (UV) renormalisation procedure for EW corrections in the 4FS are explained in Sec. 2.1. The IR singularities, in the MG5_aMC framework, [61], are dealt with via the FKS method [62,63], which is automated in the module MadFKS [64,65]. From the FKS side, the usage of the 4FS in conjunction with EW corrections does not pose additional difficulties; bottom quarks are massive and do not give rise to any collinear singularity. In practice, they are treated in the same way as top quarks. Initial-state bottom quarks are not present and IR divergences can be regulated according to the procedure explained in Sec. 3 of Ref. [52].

Complex mass scheme renormalisation and virtual matrix-element evaluation in the 4FS
In order to handle the intermediate resonances appearing in the Feynman diagrams, we adopt in our calculation the well-known complex mass scheme [66,67], in which one has to modify the renormalisation conditions yielding complex-valued renormalised parameters. Beyond LO, the carry out of the complex renormalisation procedure becomes subtle and requires very careful treatments. In particular, Ref. [52] has explored several important issues related to the computations at NLO in general. However, before this present publication, in the MG5_aMC framework NLO EW corrections and more in general complete-NLO predictions could be performed only in the 5FS, i.e., with massless bottom quarks. In order to perform the calculation of the complete-NLO predictions of any SM process in the 4FS, in particular Hbb as a case study in the present paper, we have extended the capabilities of MG5_aMC. Therefore, in this work we report also the new feature of MG5_aMC: NLO EW corrections and more in general complete-NLO predictions can also be obtained in the 4FS. Technical difficulties, as already pointed out in Ref. [52], are related to the presence of very different scales in the process and the use of the complex renormalisation conditions in it. The main concern here is that, since the complication arises from multiple-scale Feynman integrals, we have to take care of the analytic continuation from the first Riemann sheet to other sheets of two-point one-loop integrals appearing in the mass and wave-function UV renormalisation counterterms in the complex mass scheme. Instead of performing Taylor expansions in the simplified version (cf., e.g., Sec. 6.6.3 in Ref. [68]), our implementation follows a more rigorous approach, the so-called trajectory approach, first proposed in Ref. [52]. The latter does not introduce additional approximations, following the original spirit of the complex mass scheme. However, the difference between the two is in general formally beyond NLO in the SM and thus numerically insignificant [69] for the SM particle's mass spectrum. On the contrast, for a general mass spectrum in, e.g., a BSM theory, the simplified version could fail to produce the correct result, while the trajectory approach always selects the correct Riemann sheets for the multivalued complex functions. The concrete realisation of the trajectory approach in our implementation has been given in the Appendix E.2 of Ref. [52]. More specifically, the numerical routines follow Eqs. (E.44)-(E.47) in that appendix. We want to stress that although the general conceptual issues have already been extensively discussed in Ref. [52], the novel aspect in the current paper is the first complex mass scheme realisation of 4FS in MG5_aMC, as well as its validation 2 in the context of the general trajectory approach. Such an implementation will be publicly available in a future release of MG5_aMC.
The evaluation of one-loop virtual matrix elements in MG5_aMC is performed in the Mad-Loop module [61,70], by using different types of techniques for Feynman one-loop integral reduction, namely, the integrand reduction (e.g., the so-called OPP [71] and Laurent-series expansion [72] methods) or tensor integral reduction [73][74][75] approaches. MadLoop is used to automatically generate the one-loop renormalised amplitudes and to evaluate them via dynamically switching among the different one-loop integral techniques by employing the public codes CutTools [76], Ninja [77,78] and Collier [79]. An independent in-house implementation of the OpenLoops optimisation [80] enables the boosting of the fast evaluations of the virtual matrix elements. Contrary to the 5FS, in the 4FS a worry may come from the possible numerical instability occurring in one-loop evaluations due to the smallness of the bottom quark. As an example, in our calculation (i.e. complete-NLO predictions for Hbb production) we find that the self-diagnostic and recovery strategies implemented in MadLoop5 [61] are already quite effective. Namely, 99.79% phase-space points have been successfully calculated by Ninja in the double precision, and the majority of the remaining unstable phase-space points have been successfully rescued by Collier and CutTools in the double precision, while only two phase-space points (amounting to one in hundred million events) have needed the quadruple precision architecture. No event has failed to be rescued.

Topologies contributing to the Hbb final state
In order to better understand the results of our calculation, it is first of all useful to describe the various topologies of the diagrams entering each perturbative order of the complete-NLO predictions, as summarised in Tab. 1. The LO 1 originates from the "genuine" Hbb production process in the 4FS, i.e., topologies that feature a bottom-Higgs coupling, such as the one depicted in Fig. 1(a) for gg→Hbb. Also contributions from the quark-antiquark initial state are present at this order, Fig. 1(b), but their contribution is much smaller than those form the gluon-gluon initial state. The LO 3 receives contributions from "genuine" Hbb production via γγ→Hbb diagrams, but the qq→Hbb diagrams dominate at this order. Indeed, not only the γγ→Hbb process is suppressed by the photon PDF, the qq→Hbb diagrams contain an additional topology: qq→ZH production with the subsequent Z→bb decay, see the diagram illustrated in Fig. 1(f). For this reason, being the Z boson typically on shell, the LO 3 contribution is not expected to be suppressed w.r.t. the LO 1 one by a factor of order (α/α s ) 2 , as one may guess from a naive α s and α power counting. On the other hand, LO 1 and LO 3 have completely different shapes at the differential level. The LO 2 instead receives contributions only from gγ→Hbb diagrams (with a "genuine" Hbb topology) and therefore it is expected to be negligible in comparison to the LO 1 and the LO 3 .
Moving to NLO, we do not list here all the possible partonic initial states, but nonetheless we comment on the topologies appearing at any NLO i order. On the one hand, NLO 1 and NLO 2 can be viewed as the NLO QCD and NLO EW corrections to the "genuine" Hbb production,

Perturbative order
Topologies Perturbative order Topologies Table 1: Topologies entering at LO, with initial states that are explicitly specified, and at NLO. As discussed in Sec. 2.3, the terms proportional to y b y t in NLO 1 , emerging from the interference of Hbb and ggF+bb topologies, are not taken into account in our calculation.
respectively. On the other hand, NLO 3 and NLO 4 can be viewed mainly as the NLO QCD and NLO EW corrections to the Z(→bb)H production, respectively. However also new topologies enter the calculation, inducing a sensitivity to new interactions. NLO 1 receives a contribution from an additional topology: gluon fusion with an emission of a gluon splitting into a bb pair, ggF+bb, giving rise to terms of order y b y t . Similar contributions are also present in the NLO 2 , where the bb pair instead emerges from a photon or Z-boson emission from the loop. At the same order, terms proportional to y b y t can be induced also by diagrams such as the one in Fig. 1(d), which has a similar topology to the one shown in Fig. 1(a), but does not depend on y b . Similarly, diagrams like the one in Fig. 1(e) can induce a sensitivity on the HW W interaction without depending on y b . In general, EW corrections can potentially induce a sensitivity on any other SM electroweak interaction 3 and in particular, in the case of the Higgs boson, on interactions different from y b , with much larger coupling constants.
One of the most important findings of our work is that also the NLO 3 term receives a contribution from an additional new topology, namely, the VBF-like diagrams that appear at this order for the real-emission processes gq→Hbbq (see a representative diagram in Fig. 1(g)); these diagrams are the 4FS counterpart of qb→Hqb contributions to VBF in the 5FS. For this reason, NLO 3 contributions are expected to be large and to have different shapes from the LO 3 ones. Similar VBF contributions are present in the NLO 4 , where instead of a gluon, a photon is present in the initial state. Moreover, the argument of the previous paragraph regarding the sensitivity on additional EW interactions introduced by EW corrections applies also to the NLO 4 . In this case, the dominant underlying tree-level topology is the V H one, Fig. 1(f), but also contributions coming from the Hbb topology, Fig. 1(b), with the gluon substituted by a γ/Z boson, are present. In conclusion, all the perturbative orders are in principle non-negligible and exhibit different shapes at the differential level.
We remind the reader that the loop-induced ggF+bb topology leads to very large contribu-tions to the Hbb final-state. As discussed in Ref. [50] these contributions involve terms proportional to y 2 t and to y b y t , where the latter, as already mentioned before, is part of the NLO 1 and originates from the interference of the ggF+bb and "genuine" Hbb topologies. Therefore, the presence of these y t -dependent terms in addition of the very large (reducible) background due to ggF plus additional light jets, makes the extraction of y b from Hbb measurement very challenging. The fact that two additional different topologies (ZH and VBF) are present and give sizeable contributions that depend on the HZZ coupling rather than on y b , dramatically complicates the extraction of y b from Hbb production at the LHC. As we will quantify better in Sec. 3, given the smallness of the Hbb production cross section already at the inclusive level, possible extra selection cuts that reduce the dependence on y t and the HV V coupling are not only difficult to design, but also end up in killing an already rare signal.
Before moving to the next section we comment on the reliability of the usage of the pure 4FS in the context of our study. As we will better explain in Sec. 2.3, the two conditions y b = 0 and m b = 0 are inconsistent when EW corrections are taken into account. Thus, the 4FS remains the only possible choice for performing the computation of complete-NLO predictions to Hbb production. However, even considering only QCD corrections, one may argue that in such a scheme the perturbative convergence is jeopardised by the presence of large logarithms of the form α n s log k (Q/m b ) with n, k > 0, where Q is a hard scale involved in the process. For Hbb production, such potential large logarithms are only of initial-state origin and arise from the initial-state gluon splitting into a bb pair and subsequent gluon emissions from the b quarks, leading to the condition k ≤ n. In our calculation, they appear both in the "genuine" Hbb topology (twice, once for each initial-state gluon) as well as in the VBF topology (only once). In the 5FS, these logarithms would be resummed and automatically incorporated in the evolution of the bottom-quark PDF, which would be the 5FS counterpart of the initial-state gluon splitting into a bb pair that is present in the 4FS diagrams. On the other hand, as it was in general shown in Refs. [43,44], the impact of these collinear logarithms is not very large unless the process is dominated by large-x dynamics. It was also shown that the typical scale Q entering the terms α n s log k (Q/m b ) in a 4FS calculation is considerably smaller than the naively expected value, i.e., the total invariant mass of the final-state (defined excluding the bottom quarks). The same argument also suggests the usage of a rather low value for the factorisation (and renormalisation) scale, as we will do in this work and specify in Sec. 3.1. In conclusion, the combination of 4FS and 5FS predictions as done in Refs. [45][46][47][48][49] is definitely important for improving the precision of the description of the contribution from the "genuine" Hbb topology, but not compulsory for providing a sensible prediction and the 4FS can be safely used for the purpose of our work.
A different kind of potentially large logarithms (of the form log k (µ R /m b ) where µ R is the renormalisation scale) instead arise from the renormalisation of m b and y b and have a huge impact, which cannot be neglected in our work, on the predictions for Hbb. We discuss this subject in the next section.

Renormalisation of the bottom Yukawa coupling and EW corrections
An important issue in the calculation carried out here is the renormalisation condition for the parameter m b , i.e., the mass of the bottom quark. In our calculation we adopt the complex mass scheme for the massive unstable particles (the W, Z and Higgs bosons and the top quark), the G µ scheme for the electroweak interactions and the MS scheme with four active flavours for α s . The remaining renormalisation condition is the one concerning m b , which enters the calculation both via the phase-space integration, since there are two on-shell bottom quarks in the final state, and via the matrix elements, especially via the y 2 b dependence. The parameter y b is the Yukawa of the bottom quark and in the SM where v is the Higgs vacuum expectation value. While in our calculation it is natural to set an on-shell condition for m b itself, in the case of y b the typical scale involved in the interaction is much larger, being of the order of the Higgs mass m H . Performing a calculation of higher-order effects of purely QCD origin, the problem can be solved by treating y b and m b as independent parameters, with y b renormalised in the MS scheme at the renormalisation scale µ R , and m b on-shell. This strategy has shown an improved convergence of the perturbative series, automatically resumming in y b large logarithms of the form log(µ R /m b ), and a reduction of the difference with the corresponding calculation in the 5FS. Moreover, it allows to independently vary y b and m b .
The problem arises from the fact that when EW corrections are calculated the SM relation y b = √ 2m b /v must be enforced 4 in order to have the cancellation of UV divergences. This fact has two consequences for a SM calculation at NLO EW or complete-NLO accuracy for Hbb production: and therefore the computation is not feasible, • in the 4FS the renormalisation condition for m b fixes the renormalisation condition for y b and vice versa.
Although the main phenomenological results of this paper are not strictly based on precision physics, namely percent or higher-level accuracy, a correct treatment of the input parameters for the LO 1 and NLO 1 is crucial in order to obtain sensible result; using /v can lead to results larger by a factor of 2 to 3 than in the case of In order to amend this situation, we adopt the following procedure. First of all, we perform the purely QCD calculation (LO 1 and NLO 1 ) by employing the MS scheme for the Yukawa of the bottom (y respectively, are consistent. Following the recommendation from Sec. IV.2.2.a of Ref. [82], we adopt a one-loop QCD transition between the two schemes. 5 At the same time, we exclude the contribution from interference terms between ggF+bb and Hbb topologies in NLO 1 , which are proportional to y b y t and have already been studied in Ref. [50] at higher accuracy and shown to have a mild impact on the cross section. Second, we combine the calculation of the LO 1 and NLO 1 contributions in this scheme, dubbed LO MS 1 and NLO MS 1 | yt=0 respectively, with the remaining complete-NLO terms where the Yukawa of the bottom is renormalised on-shell /v. In doing so we do not simply add the different perturbative orders, but we combine them in a sort of multiplicative approach, namely, we first define and then All quantities without a superscript in Eqs. (3) and (4) are meant with m b renormalised on shell, at variance with those with MS . 4 This fact is strictly true in the SM. For example, in an effective-field-theory approach the relation y b = √ 2m b /v can be modified while keeping the possibility of performing EW corrections, e.g., by adding a dimension-6 operator of the form Φ † Φ Λ 2QL Φb R , where Φ is the Higgs doublet. However, the renormalisation conditions for this kind of calculation is much more involved, see e.g. Ref. [81]. 5 The reason is that renormalons are present in the quantity and the perturbative series does not converge, as discussed in Ref. [48]. Therefore, at variance with a general convergent series, the common lore "the higher is the precision the better it is" does not apply here.
Beside the renormalisation of the Hbb vertex, the NLO 2 term contains several other types of contributions, such as EW Sudakov logarithms that depend on y b only via the underlying LO 1 contribution and therefore can be naturally rescaled by a factor (m MS b /m pole b ) 2 in order to provide an improved prediction; the term NLO MS 2 in Eq. (3) precisely corresponds to the rescaling of the NLO 2 contribution by this factor. In principle one may think also to add a term NLO 2 × (NLO MS 1 | yt=0 /LO 1 ) as typically done in the multiplicative combination of NLO QCD and EW corrections. We have checked this alternative approach and found negligible differences with the results obtained via Eq. (4). Finally, we remind the reader that the orders LO 3 , NLO 3 and NLO 4 are dominated by ZH and VBF configurations, for which the issue concerning the renormalisation of the Hbb interaction is not relevant.
For simplicity, in the rest of the paper we will use also the notations defined in the following 3 Numerical results

Input parameters
We provide numerical results for proton-proton collisions at the LHC with a center-of-massenergy of 14 TeV, as planned for the Run-III and the High-Luminosity run [83]. We perform the calculation using the complex mass scheme and the following input parameters 6 where we have set Γ H = 0 since there is an external on-shell Higgs in our calculation. EW interactions are renormalised in the G µ -scheme with We set the pole mass of the bottom quark to when the difference between the two schemes is evaluated at one-loop level, as motivated in Sec. 2.3. We do not evaluate uncertainties related to the value of m b ; they are discussed, e.g., in Refs. [48,49] and are of the order of a few percents.
For the factorisation scale µ F and the renormalisation scale µ R , which enters the definition of α s (µ R ) and y b = √ 2m MS b (µ R )/v in the NLO QCD calculation, we use a central value where the index i runs over all the final-state particles. The µ R dependence of α s is directly taken from the PDF set employed in the calculation, while we evolve m MS b (µ R ) at four loop in QCD [84,85]. Scale uncertainties are evaluated by independently varying the factorisation and renormalisation scale in the range µ 0 /2 < µ F , µ R < 2µ 0 .
Phase-space integration is performed with no constraints on the b-quark momenta. On the other hand, we will provide results for the full phase space as well by setting constraints on the number of b-jets and possibly light jets. Jets are clustered with the anti-k T algorithm [86] as implemented in FastJet [87], with the distance parameter R = 0.4. Jets are required to have p T > 30 GeV and pseudorapidity |η| < 4.5. Photons are clustered into jets and therefore also a real emission of a single photon can form a separate jet. 7 In the case of jets that have at least a bottom quark or antiquark among their constituents, the requirement |η| < 2.5 has to be satisfied in order to be classified as b-jets, otherwise they are tagged as light jets. As we will discuss in Sec. 3.2, we will also explore the effects of a jet-veto in the entire phase-space region |η| < 4.5.
Since we calculate NLO EW corrections, as discussed in Sec. 2.2, processes featuring initialstate photons are present and therefore a parameterisation of the photon PDF is necessary. Furthermore NLO QED effects in PDFs evolution have to be taken into account. Even more important, the PDFs should be in the 4FS in order to correctly take into account QCD effects. However, to the best of our knowledge, at the moment there are no public PDF sets including NLO QED effects and a photon density in the 4FS. For this reason, we have checked the numerical impact of the photon PDF and NLO QED effects by comparing results obtained via the usage of the PDFs set NNPDF31_nnlo_as_0118_luxqed [88], which is based on the PDF fit NNPDF3.1 [89] and the photon parameterisation of LUXqed [90,91], and the PDF set NNPDF31_nnlo_as_0118, which includes neither a photon PDF nor NLO QED effects in the evolution. These two PDF sets are both in the 5FS, but this is irrelevant for the sake of estimating QED effects in the PDFs. We find that effects related to the QED evolution are of the order of 1% for LO 1 , while the LO 2 term, which is purely γg-induced, gives a contribution that is 0.1% of the LO 1 one. Therefore, neglecting QED effects in PDFs is completely justified for the study carried out in this work. We conclude that we can safely use the PDF set NNPDF31_nnlo_as_0118_nf_4l which does not account for QED effects, but is designed for calculations in the 4FS. 8 Nevertheless, we advocate the necessity to have public PDF sets including QED effects in the 4FS.

Inclusive results
We now turn to the presentation of results, starting from total cross sections for Hbb production at 14 TeV defined for different jet categories. In Tab. 2, we list predictions computed at different perturbative accuracies, according to the definitions in Eqs. (5)-(9). We show results for different selection cuts on b-jet multiplicities, namely, • NO CUT: No restriction on the momenta of the final-state particles, • N j b = 1: Exactly one b-jet, with and without a veto on light jets, 7 We remind the reader that in many LHC analyses a jet is defined with up to 99% of its energy of electromagnetic origin and up to 90% of it that can be carried by a single photon. See Ref. [56] for further details on this subject. 8 In principle, one could use the PDF set NNPDF31_nnlo_as_0118_luxqed and remove the impact of the fifth flavour from the running of α s and to the DGLAP equation for the PDF evolution in order to be consistent at NLO QCD accuracy. For instance, one may adopt a strategy similar to the one explained in Sec. 2.2 of Ref. [92], which is in turn based on Ref. [93]. However, logarithms of the form log(µ F,R /m pole b ) would be present and not resummed, especially when varying the scale.
• N j b ≥ 1: At least one b-jet, with and without a veto on light jets, At our accuracy, complete-NLO, there cannot be more than two b-jets and therefore N j b ≥ 2 ⇐⇒ N j b = 2. Numbers in parentheses refer to the case where the light-jet veto is applied. In the second column we show total rates for the central scales µ F = µ R = µ 0 , together with relative scale uncertainties, while in the third column we show the corresponding ratios with the central-scale LO QCD predictions. The relative impact of the different perturbative orders is further documented in Tab. 3, where the ratios of all the different contributions entering the complete-NLO predictions (NLO all ) divided by LO QCD are separately displayed, see also Eqs. (3),(4) and (9). We recall that LO 2 is exactly zero since we use a PDF set without a photon density, and therefore its contribution is not displayed in Tab. 3. We also remind the reader that the term LO 1 is equivalent to LO QCD , but with the Yukawa of the bottom renormalised on-shell,

Description of the results
We start the discussion of the numerical results by commenting the numbers in Tabs. 2 and 3 obtained without applying the light-jet veto. We will then move to the case with the light-jet veto and finally we will draw our phenomenological conclusions in Sec. 3.2.2: at variance with the naive expectation, the measurement of total rates for Hbb production is not leading to a direct sensitivity to y b , regardless of the selection cuts that are used.
Results without the light-jet veto The most important feature that can can be observed in Tab. 3 is that the relative impact of LO 3 , NLO 3 and NLO 4 grows with N j b . First of all, these contributions, before setting cuts, are not dominated by the "genuine" Hbb topology, but rather by the ZH and (except LO 3 ) VBF topologies. Then, while in the "genuine" Hbb topology with the gluon-gluon initial state ( Fig. 1(a)), which dominates LO QCD , NLO MS 1 and NLO MS 2 , both the bottom quarks tend to be collinear to the beam-pipe axis, in the VBF topology this holds true for only one of the two bottom quarks and for none of them in V H. Therefore, the probability that a bottom quark b is either soft or falls outside the rapidity region in which b-jets are tagged, |η(j b )| < 2.5, is higher for the "genuine" Hbb topology than for ZH and VBF topologies. The same behaviour has been observed in Ref. [50] regarding the comparison with the ggF+bb topology. The net effects is the aforementioned growth of the relative impact of LO 3 , NLO 3 and NLO 4 with N j b .
The same feature can be observed also in Tab. 2 by comparing the LO QCD , NLO QCD and NLO QCD+EW predictions, which do not include the LO 3 , NLO 3 and NLO 4 contributions, with the LO and NLO all ones, which do include (part of) them. In fact, according to Eqs. (5)-(9), since we set the photon PDF to zero, we exactly have LO = LO QCD + LO 3 and NLO all = NLO QCD+EW +LO 3 +NLO 3 +NLO 4 . Therefore, as already demonstrated in Refs. [58,60,94,95] for other processes, contributions formally suppressed by the (α/α s ) naive power counting can actually be numerically much larger than expected, especially when specific phase-space cuts are imposed. We remind the reader that each of the rates for N j b ≥ 1 in Tab. 2 is equal to the sum of the corresponding ones for N j b = 1 and N j b ≥ 2. By looking at the numbers for N j b ≥ 2 one can understand the large difference between the case N j b = 1 and N j b ≥ 1. With N j b ≥ 2 the complete-NLO prediction, NLO all , is 10.8 times larger than the LO QCD one. The LO 3 is 7.9 times larger than the LO QCD , the NLO 3 is 2.4 times larger, and the NLO 4 is -60% of the LO QCD . As an example, via a naive (α/α s ) power counting the NLO 4 would be expected to be of the order of 0.01% of the LO QCD . Although smaller in size, a similar pattern is observed also for the case N j b = 1 and therefore also for the case N j b ≥ 1. One can also notice that moving    from N j b ≥ 1 to N j b = 1, the LO 3 contribution is strongly reduced, roughly by a factor of 11, while the NLO 3 is reduced much less, roughly by a factor of 2.5. This is a clear sign that the contribution of the VBF topology to the NLO 3 is sizeable. While the ZH topology tend to have two separate b-jets, the VBF one mostly exhibits a bottom-quark collinear to the beam-pipe axis and the other one sufficiently central in order to form a b-jet. Therefore, once the N j b ≥ 2 contribution is removed, only the ZH topology is strongly suppressed. This argument will be corroborated by the analysis of the m(j b,1 , j b,2 ) distribution, i.e., the invariant mass of the two b-jets, which is presented in Sec. 3.3.
Regarding the NLO MS 2 contribution, i.e. what is typically denoted as the NLO EW corrections, it is of the size expected by the naive (α/α s ) power counting: of the order of a few percents of LO QCD predictions. Moreover, it mildly depends on the value of N j b . The reason is that at this order there are no new topologies opening, at variance with the NLO 3 and NLO 4 cases. If we had not consider the quantity NLO MS 2 , as defined in Eq. (3) (see also Eq. (5)), but directly NLO 2 from Eq. (2), the contribution of NLO EW corrections would have been larger. Indeed, as can be seen in Tab. 3, the ratio (LO QCD /LO 1 ) is ∼2.4. This ratio has a small dependence on N j b that is induced by the renormalisation scale of y b , which is dynamical (see Eq. (15)) and therefore induces not only a global rescaling w.r.t. the LO 1 term, which has been calculated with on-shell y b , but also mild differences in shapes. As already mentioned, NLO EW corrections have already been calculated in Ref. [51]. However, at variance with Ref. [51], not only we identify NLO EW corrections as the NLO MS 2 term rather than simply the NLO 2 one, but we also include all the possible initial states contributing to this order. In Ref. [51], only the gluon-gluon initial state has been considered.
NLO QCD corrections, namely the NLO MS 1 | yt=0 term, have already been calculated in the past [24,25] and are sizeable. Still, with the exception of the case "NO CUT", they are in general smaller than the NLO 3 and LO 3 contributions. On the other hand, the NLO MS 1 | yt=0 term is especially relevant for what concerns scale uncertainties. While the LO QCD predictions have relative scale uncertainties of the order ∼ +50% −30% , NLO QCD predictions have relative scale uncertainties of the order 15-20% and even smaller for the N j b ≥ 2 case. If we had not implemented the MS scheme for y b , scale uncertainties would had been smaller at LO in QCD (LO 1 ), since y b would not depend on µ R , and also at NLO in QCD (LO 1 + NLO 1 | yt=0 ). However, this reduction of scale uncertainties should be interpreted as an underestimate of higher-order effects by the use of y b in the on-shell scheme rather than a more accurate prediction. Concerning the NLO MS 2 term, its impact on scale uncertainties is below the 1% level, as it can be seen by comparing NLO QCD and NLO QCD+EW predictions. Instead, moving from NLO QCD+EW to NLO all predictions, the size of the scale-uncertainty band decreases in any N j b category. The reason is that the LO 3 contribution has a much smaller scale dependence w.r.t. the LO QCD one, since at this order the ZH topology does not depend neither on y b nor on α s ; its scale dependence originates only from PDFs and thus from µ F . This can be seen by comparing the LO QCD predictions with the LO ones, where the latter are exactly equal to the former plus the LO 3 contribution. The NLO 3 contribution introduces a µ R dependence via the presence of one power of α s , but it also further reduces the dependence on µ F . Altogether, these effects lead to the reduction of the size of the scale-uncertainty band from NLO QCD+EW to NLO all .
Results including the light-jet veto We now comment the results where the veto on light jets is applied, namely, the number of Tabs. 2 and 3 that are in parentheses. First of all, it is worth to notice that the light-jet veto affects also LO i results because b-jets are tagged only in the |η(j b )| < 2.5 region. When 2.5 < |η(j b )| < 4.5 the jet is actually tagged as light and therefore the light-jet veto has an effect on it. Moving to NLO QCD , NLO QCD+EW , and NLO all predictions, the first comment about them is that scale uncertainties for results with the jet veto do not largely increase w.r.t. the corresponding cases without it, rather they mildly decrease. This is a clear sign that jet-veto resummation or the matching with the shower effects is not mandatory for obtaining sensible results. The situation is slightly different in the case N j b ≥ 2, where we have observed much larger scale uncertainties and therefore we have omitted them in Tabs. 2 and 3. The case "NO CUT", without the light-jet veto, has been reported in Tabs. 2 and 3 in order to document the result of our calculation and better interpret the N j b categorisation. On the other hand, we already know that its contribution is about 100 times smaller than inclusive ggF production and therefore not suitable for a sensitivity-study on Hbb production and especially on y b . For this reason, we have chosen to not show the case of a light-jet veto and N j b ≥ 0, and in conclusion we consider the light-jet veto option only for the cases N j b = 1 and N j b ≥ 1.
Like in any fixed-order calculation, the light-jet veto has a sizeable impact on the NLO QCD K-factor, i.e., the σ NLO QCD /σ LO QCD ratio, as can be seen in Tab. 2. Non-negligible effects are present also for the LO 3 and therefore the LO predictions, as can be respectively seen in Tabs. 3 and 2. However, the largest impact of the light-jet veto is on the NLO 3 contribution and therefore the NLO all predictions. While without the light-jet veto the NLO 3 contribution is of the same size of the LO QCD one for both the N j b ≥ 1 and N j b = 1 categories, applying the light-jet veto the (central value of the) NLO 3 contribution almost vanishes in the case of N j b ≥ 1 and drops to only ∼ 20% of the LO QCD one when N j b = 1. The reason is that the VBF topology typically has one light-jet induced by the light quark in the final state and possibly one additional light-jet due to one of the two bottom quarks, which is usually at large rapidities. Therefore the veto has a huge effect on the contribution from this topology. Moreover, the NLO 3 has a large contribution from "QCD corrections" to the ZH topology, which includes gluon emissions from the bottom quarks from the Z decays. The light-jet veto has a large impact also on these configurations, especially in the case of N j b = 2, which is present also in N j b ≥ 1. This is the reason why the effect of the light-jet veto on the NLO 3 contributions is slightly larger in the case N j b ≥ 1 than in the case N j b = 1. As a last remark, we notice that the impact of the light-jet veto is instead negligible on NLO MS 2 and NLO 4 contributions.

Prospects on the y b measurement
On the basis of the previous discussion and of the results of Tabs. 2 and 3, we now comment on what are the prospects of a direct determination of y b via the Hbb measurement at the LHC. For the sake of clarity, in the following discussion we will associate specific perturbative orders to specific Higgs couplings: where adopting the κ-framework notation [96] we denote the HZZ interaction as κ Z . Relations (16)-(21) also imply Clearly, as also pointed out in Sec. 2.2, the NLO MS 2 and NLO 4 terms involve contributions that depend on additional couplings and that can even not depend at all on y b and κ Z , respec-  Table 4: Fraction of the cross section scaling as y 2 b for different phase-space cuts. The first column is based on the results from our calculation in Tab. 2. The second column is based on results from Ref. [50]. The third column is based on the numbers in the first and second column. Details are explained in the text.
tively. However, one can understand from the discussion of Sec. 3.2.1 that the numerical impact of NLO MS 2 and NLO 4 terms, and therefore of such contributions, is negligible w.r.t. the other perturbative orders involved in the calculation. Moreover, as it will become more clear in the following, taking into account a more realistic and more complex coupling structure in a given perturbative order would make our argument even stronger. In other words, relations (16)- (24) are devised for simplifying the discussion, but our conclusions do not depend on them.
For the same N j b categories of Tabs. 2 and 3, in the first column of Tab. 4 we report the ratio of the NLO QCD+EW and NLO all predictions, here denoted as σ NLO QCD+EW and σ NLO all . Both of them are our best predictions for respectively the O(y 2 b ) cross section, denoted in the following also as σ(y 2 b ), and the sum of it with the O(κ 2 Z ) cross section, denoted in the following also as σ(κ 2 Z ). Via the ratio σ NLO QCD+EW /σ NLO all we can determine the fraction of the measured cross section that actually depends on y b . Once again, we remind the reader that the case "NO CUT" is purely academic, since the signal from inclusive ggF Higgs production exceeds the one of Hbb production by a factor of 100. Thus, one needs to tag at least one b-jet and we already know that also after that the ggF+bb contribution is large, so we should at least suppress the ZH and VBF topologies, which yield σ(κ 2 Z ). The category N j b ≥ 2 has very small rates (see Tab. 2) and the lowest σ NLO QCD+EW /σ NLO all ratio, due to the large contribution of the ZH topology, therefore it is not expected to be the best option in order to gain sensitivity on y b . This also explains why the category N j b = 1, which does not include N j b ≥ 2, has a larger σ NLO QCD+EW /σ NLO all ratio w.r.t. the category N j b ≥ 1, which does include it. However, in both the N j b = 1 and N j b ≥ 1 categories, the VBF contribution is still large, but the light jet-veto (numbers in parentheses) helps in reducing it. In conclusion, the best option seems to be the N j b = 1 category with a light-jet veto, where 60% of the signal depends on y b .
So far, however, we have completely neglected the contribution of the ggF+bb topology, which leads to O(y b y t ) contributions, σ(y b y t ), and especially O(y 2 t ) contributions, σ(y 2 t ). In order to amend this situation we use the results of Ref. [50], where σ(y 2 b ), σ(y b y t ), and σ(y 2 t ) have been calculated at NLO QCD accuracy. Using the numbers of Tab. 1 in Ref. [50], in the second column of Tab. 4 we report the ratio of the cross section calculated including only the "genuine" Hbb topologies or adding also the ggF+bb one. In other words, σ(y 2 b ) divided by σ(y 2 b )+σ(y b y t )+σ(y 2 t ). As can be seen, the impact of σ(y b y t ) and σ(y 2 t ) is huge and therefore cannot be neglected for our purposes.
The same definitions of b-jets have been used in Ref. [50] and in our work. A few differences in the input parameters are present, but their impact is expect to be minor, especially when ratios of cross sections are considered. In particular, we have explicitly checked that the difference for the collision energy, 13 TeV in Ref. [50] and 14 TeV in the present work, has little effect on the ratios. The only results that we cannot derive from Ref. [50] are those for the case with a light-jet veto. On the other hand, for the case without a light-jet veto, we can combine the results from the first and second column of Tab. 4. Since in the second column we have σ(y 2 b ) divided by σ(y 2 b )+σ(y b y t )+σ(y 2 t ), if we assume that the first column is σ(y 2 b ) divided by σ(y 2 b )+σ(κ 2 Z ), we can derive σ(y 2 b ) divided by σ(y 2 b )+σ(κ 2 Z )+σ(y b y t )+σ(y 2 t ), which is the quantity displayed in the third column. The result is striking: in none of the realistic N j b categories σ(y 2 b ), i.e., the component of the total cross section that scales as y 2 b , is larger than 16%. As we will see in the next section, differential information is also in general not helping in improving this ratio. Also, the light-jet veto option cannot substantially alter this picture, as can be seen by the number in the first column of Tab. 4, so this option can also be safely ruled out.
We want to stress that, if we consider σ(y 2 b ) as the "signal" in an experimental analysis, in this work we are not considering a realistic comparison between the signal and its backgrounds. At this stage, regarding the backgrounds, we are considering only the irreducible backgrounds, without even taking into account the Higgs boson decays. Needless to say, if we had taken into account also the irreducible and reducible backgrounds for a given signature that is induced by a specific Higgs-boson decay, the situation could have only got worse. From the theoretical side, the same applies if instead of assuming the simplified relations (16)-(24) we would have taken into account the complete coupling dependence. In the next section, we will explore the last hopes of identifying phase-space regions where the sensitivity on σ(y 2 b ) may be strongly enhanced. We can anticipate, that this is not the case.

Differential distributions
We start the discussion about differential distributions by analysing the m(j b,1 , j b,2 ) observable, which can be obtained in our analysis only for N j b ≥ 1 and N j b ≥ 2 and is exactly the same in the two cases, since m(j b,1 , j b,2 ) is defined only for N j b = 2. By looking at this distribution we can definitely prove that the NLO 3 order is populated by large contributions from the VBF topology, beside the ZH one. After that, we will consider many more observables for the cases N j b ≥ 1 and N j b = 1.
In Fig. 2, we show the m(j b,1 , j b,2 ) distribution at different accuracies (LO QCD , LO, NLO QCD , NLO QCD+EW , NLO all ) together with their scale uncertainties. The left plot refers to the case where the light-jet veto has not been applied, while in the right one we show results with the light-jet veto. In each plot, we show in the lower inset the same quantities of the main panel normalised to the central value of the NLO QCD prediction.
As can be seen in Fig. 2, the m(j b,1 , j b,2 ) ∼ m Z region is completely dominated by the LO prediction, which contains the LO 3 contribution, the one involving the ZH topology. The NLO 3 contribution, which is contained in the NLO all prediction, involves QCD corrections to the ZH topology, such as the emission of gluons from the bb pair stemming from the Z boson decay. The radiation of gluons form the b quarks together with the presence of the Z resonance leads to a large amount of events migrating from the m(j b,1 , j b,2 ) ∼ m Z region to smaller values of m(j b,1 , j b,2 ). This behaviour is typical for any invariant mass distribution of decay products of a resonance, when either QCD or QED emissions are considered. However, at variance with this standard picture, in the left plot of Fig. 2 we can see that the difference between the NLO all and LO prediction, which is mainly induced by the NLO 3 contribution, is large also for m(j b,1 , j b,2 ) m Z . This effect is precisely induced by the presence of VBF configurations, which on the other hand are suppressed when a light-jet veto is applied, as can be seen in the right plot. In Tabs. 2 and 3, we did not show results with the light-jet veto for N j b ≥ 2 since scale uncertainties are too large. Indeed, this feature can be seen in the right plot. The analysis of the m(j b,1 , j b,2 ) spectrum shows also that even applying a cut around the m(j b,1 , j b,2 ) = m Z value, the result would be still contaminated by VBF configurations.
In the right plot the light-jet veto is applied.
We now proceed to the analysis of differential distributions for several observables in the case N j b ≥ 1 and N j b = 1. First of all, beside documenting the result obtained, we want to explore the possibilities of enhancing the sensitivity on the σ(y 2 b ) contribution. In each of the Figs. 3-8 we show distributions for a specific observable for N j b ≥ 1 (upper plots) and N j b = 1 (lower plots) without (left plots) and with (right plots) the light-jet veto applied. We consider the following distributions: the transverse momenta and the pseudorapidity of the hardest bjet, p T (j b,1 ) and η(j b,1 ), the transverse momenta and the rapidity of the Higgs boson, p T (H) and y(H), the absolute value of the difference of the Higgs and hardest b-jet pseudorapidities, |∆η(H, j b,1 )|, and finally the transverse momenta and the pseudorapidity of the light-jet, p T (j l ) and η(j l ). Since the last two observables are not defined in the case of the light-jet veto, we combine them in Fig. 8.
All the plots of Figs. 3-8 have the same layout of those in Fig. 2, which has already been described. First of all, one can see that also at the differential level the NLO MS 2 contribution, which is equal to the difference between the NLO QCD+EW and NLO QCD predictions, is negligible. In absolute value, it reaches at maximum few percents of the NLO QCD prediction in the tails of the transverse-momentum distributions. For this reason, the NLO all /NLO QCD ratio can be interpreted as the differential version of the ratio [σ(y 2 b ) + σ(κ 2 Z )]/σ(y 2 b ), namely the inverse of the quantity displayed in the first column of Tab. 4. The higher is the NLO all /NLO QCD ratio, the smaller is the fraction of the cross section that depends on y b . The most important result that can be obtained by the analysis of all these plots is that whenever we look at phase-space regions that do not correspond to the bulk of the cross-section (large values of p T , |η| or |y|, etc.), the NLO all /NLO QCD ratio increases. In other words, applying cuts that depend on any of the observable we have considered, total rates diminish and at the same time the fraction of the cross section that depend on y b decreases. The only exception is the |∆η(H, j b,1 )| distribution, especially when the light-jet veto is applied. However, in order to halve the σ(κ 2 Z ) term and bring it to roughly 30-40% of σ(y 2 b ), rates have to be suppressed by a factor of 10. Thus, no real improvement can be gained. In conclusion, the sensitivity on y b cannot be improved even via the information at the differential level.
Although the main message of our phenomenological analysis has already been conveyed, we now report the other important features of plots in Figs. 3-8. We start with the p T (j b,1 ) distribution in Fig. 3. By comparing NLO QCD and LO QCD predictions one can see that the relative impact of the NLO MS 1 | yt=0 contribution is rather flat if the light-jet veto is not applied, both in the N j b ≥ 1 and N j b = 1 cases. By applying the light-jet veto, the NLO MS 1 | yt=0 term becomes negative at large p T (j b,1 ) values, with a larger impact for the case N j b ≥ 1. Both with and without the light-jet veto, NLO QCD scale uncertainties are much smaller than the LO QCD ones, also at the differential level. As already said, the NLO MS especially for N j b ≥ 1 and/or applying the light-jet veto. It is important to note that in the case of the light-jet veto this effect is due to the suppression of the NLO QCD prediction; the LO prediction is only mildly affected by the light-jet veto also for large p T (j b,1 ). Moving to the NLO all prediction, which in particular includes the NLO 3 term, also this quantity is larger than the NLO QCD prediction, and also the NLO all /NLO QCD ratio grows for large p T (j b,1 ) values, especially for N j b ≥ 1. On the other hand, the impact of the light-jet veto is the opposite than in the LO case; the NLO all prediction is strongly reduced, especially at large p T (j b,1 ) values. This is not surprising, since the VBF topology typically displays a light jet and therefore is completely suppressed.
In the case of the η(j b,1 ) distribution, Fig. 4, similar considerations to the ones discussed for the p T (j b,1 ) distribution apply. The only difference is that the LO QCD /NLO QCD ratio is flat,  Fig. 3 apply also here. The situation is instead different for p T (H) ≤ 30 GeV. Indeed, the LO QCD prediction is smaller than the NLO QCD one and strongly decreases close to the threshold, especially for the case N j b = 1. This is a pathological behaviour that is typical of fixed-order calculations in the presence of hard cuts. 9  Besides, the condition p T (b 2 ) < 30 GeV and/or |η(b 2 )| > 2.5 must be satisfied, otherwise b 2 would form another b-jet. These requirements all together pose strong constraints on the b 2 phase-space, especially for p T (H)→0, suppressing the LO QCD and LO predictions. By adding a new particle in the final state, as in any NLO prediction, these hard cuts are removed and the pathological behaviour disappears. We also notice that the LO prediction, not the LO QCD one, considerably increases in this region moving from the N j b = 1 to N j b ≥ 1 case. This is due to the presence of the ZH topology in the LO 3 term. By allowing more than one b-jet, the LO 3 logical behaviour for p T (H) ≤ 30 GeV.  Figure 7: The |∆η(H, j b,1 )| distribution for N j b ≥ 1 (up) and N j b = 1 (down). In the right plots the light-jet veto is applied. and in turn LO predictions can easily satisfy the relation p T (j b,1 ) + p T (H) + p T (b 2 ) = 0. Indeed, the b quarks emerging from the Z boson decay are back-to-back in the Z boson rest-frame and not so rarely with p T (b) > 30 GeV. This leads to the presence of two b-jets, N j b = 2, which does not suppress so much the LO 3 contribution and in turn the LO contribution w.r.t. the N j b = 1 case, as can also be seen in Tab. 2. Instead, in the case of LO QCD , which is dominated by "genuine" Hbb topology, bottom quarks are typically emitted collinearly to the beam-pipe axis.
In principle, also for the LO QCD case, the conditions p T (H) < 30 GeV, p T (j b,1 ) > 30 GeV and N j b ≥ 1 could be satisfied when N j b = 2. In practice, at variance with the LO case, at LO QCD this condition leads to large suppressions of the cross sections, as can also be seen in Tab  than in the NLO QCD or LO case. On the other hand, we notice also that while the NLO QCD and the NLO all predictions are reduced by the light-jet veto, again the LO one is not.
In the case of the y(H) distribution, Fig. 6, the most important feature is the growth of the LO/NLO QCD and NLO all /NLO QCD ratios for large |y(H)| values, especially if the light-jet veto is not applied. This can be understood by the fact that the Higgs boson recoils against the bb pair and possibly an additional real emission. Therefore, at the partonic level, i.e., before the convolution with the PDFs, at LO or LO QCD accuracy and for large values of y(H) we have |y(H)| ∼ |η(H)| = |η(bb)| ∼ |y(bb)|. However, while the bb pair stems from the Z boson decays in the ZH topology, and therefore the entire bb pair tends to move in the same direction, in the case of a boosted Z boson, in the "genuine" Hbb topology the bottom quarks tend to be emitted collinearly to the beam-pipe axis and back-to-back to each other. Therefore, in the LO predictions, and especially in the NLO all one which can get a further boost from the real emissions, the large y(H) region is more populated than in the NLO QCD predictions. The light-jet veto reduces this effect, clearly more for the NLO all case. Figure 7 shows the |∆η(H, j b,1 )| distribution, which as we have already said is the only one that displays a reduction of the NLO all /NLO QCD (and also LO/NLO QCD ) ratio by moving away from the bulk of the cross section, i.e., going towards large |∆η(H, j b,1 )| values. For small |η(H)| values, where the cross section is the largest, the probability of having the hardest b-jet with small |η(j b,1 )| values is higher for the ZH topology (the LO 3 and LO contributions) than in the "genuine" Hbb topology (the LO QCD contributions), since in the latter bottom quarks tend to be emitted collinearly to the beam-pipe axis and back-to-back to each other. Also, at large |η(H)| values, in the case of the ZH topology in the LO 3 the Higgs boson mostly recoils against the bb pair with the bottom quarks moving in the same direction, while in the "genuine" Hbb topology of LO QCD the two bottoms tend mostly to have opposite directions, leading to one of them having (in the partonic rest frame) the pseudorapidity larger than the one of the Higgs boson in absolute value and with opposite sign. This dynamics is the origin of a flatter |∆η(H, j b,1 )| distribution for LO QCD and NLO QCD predictions w.r.t. the LO ones. The presence of real emissions and the VBF topology flattens the distribution moving from LO to NLO all accuracy. The flattening is even stronger moving from the N j b ≥ 1 to the N j b = 1 case, which reduces the ZH contribution.
We finally discuss Fig. 8, which displays the p T (j l ) distribution on the left and the η(j l ) one on the right. In the case of p T (j l ), we see how going to large p T (j l ) values, the LO contribution decreases w.r.t the NLO QCD one. We recall the at LO the light jets are only given by bottom quarks with pseudorapidity larger than 2.5 in absolute value, while in the NLO QCD predictions they can be genuine light-jets, with no b-quark inside them. Therefore, by requiring large pseudorapidities it is more difficult to achieve large transverse momenta. On the contrary, going to large p T (j l ) values, the NLO all contribution increases w.r.t NLO QCD one. Indeed, the light-jet in the VBF topology would not diverge in the limit p T (j l )→0, at variance with those from real QCD (or QED) emissions. For this reason the p T (j l ) spectrum at NLO all is much flatter than the one at NLO QCD accuracy. Moving to the η(j l ) distributions, the right plots clearly display the fixed-order pathological behaviour for this observable in our calculation. Indeed, in the region |η(j l )| < 2.5, the LO QCD and LO predictions are exactly equal to zero. This is the reason why NLO QCD and NLO all scale uncertainties are smaller outside of this region; in the range |η(j l )| < 2.5 are in fact "LO-type" predictions. It is interesting to note how the peak of the distribution at LO and NLO all accuracy is in the region |η(j l )| 2.5, so the pseudorapidity coverage of the b-jet tagging has a non-trivial impact on the numbers obtained in our work.

Conclusions
The precise measurement of the Higgs boson couplings is one of the major goals of the LHC program. In particular, in the case of the Higgs-fermion Yukawa couplings, this translates in the need of measuring the relevant production mechanisms and/or decay modes of the Higgs boson. For what concerns the bottom quark, in principle the extraction of y b from the measurement of Hbb production would be subject to less theoretical assumptions than the corresponding H → bb decay, whose branching ratio depends on all the other Higgs decay modes. However, the measurement of Hbb production is plagued by various backgrounds, with very large rates.
In this paper, by computing the QCD and EW complete-NLO predictions for Hbb production in the 4FS, we have shown that the irreducible backgrounds involving the HZZ interactions completely submerge the "genuine" y 2 b -dependent Hbb signal. Among these backgrounds, one has both contributions where the accompanying bottom quarks originate from a resonant decay, ZH production with Z → bb, but also contributions with a non-resonant spectrum, namely the b-initiated VBF topology. Both these classes of irreducible backgrounds have very large cross sections when compared to the "genuine" Hbb signal, and because of the different underlying structures, it is extremely complicated to reduce their rate without de facto killing also the signal. On this respect, in our study, we have considered different set of cuts, both on the number of b-tagged jets and possibly vetoing extra light-jet radiation; in all cases the aforementioned backgrounds are at least as large as the signal.
Once the other irreducible backgrounds are also taken into account, namely those coming from the ggF+bb topology depending on y t , their sum overwhelms the signal by about one order of magnitude. Thus, we find that it is tremendously difficult, if not impossible, to directly extract the bottom-quark Yukawa y b via the measurement of SM Hbb production at the LHC. Unless y b is significantly enlarged by new physics, a scenario which is strongly disfavoured by the H→bb experimental measurements, even for BSM scenarios the direct determination of y b via this process seems to be hopeless at the LHC. We have also investigated several differential distributions and we have found that moving away from the bulk of the cross section, not only the rates but also the fraction of them that depends on y b decreases.
We reckon that our study is performed at fixed-order, and neglects parton-shower, hadronisation and detector effects. Taking into account these effects it would be possible to perform a more realistic simulation. However, doing so, one should also consider on top of the Hbb irreducible backgrounds also those for the targeted Higgs decays and especially the reducible backgrounds. In general, this will further reduce the signal-over-background ratio. Thus, we do not expect our conclusions to be altered, rather reinforced. Beside the case of the measurement of "genuine" Hbb production, results presented in this paper can also be relevant for the background estimation for other Higgs processes, in particular HH production with one Higgs boson decaying to b quarks. Similarly, the Hbb final state has to be taken into account when precise predictions for inclusive Higgs boson production are calculated. Therefore, the calculation presented in this paper can be in principle exploited also for this purpose, however, in this case one should pay attention to not double-count ZH and VBF contributions, which can be computed at higher accuracy via dedicated calculations. We leave studies in this direction for future work. Moreover, regardless of our phenomenological findings, the Hbb process remains a key process for the improvement and understanding of techniques for the computation of higher-order corrections in QCD, theoretical developments for the combination 4FS and 5FS computations at different perturbative orders and, as shown in this paper, for a better insight in the structure of the renormalisation condition in the EW sector.
Finally, beside the phenomenological results, we have extended the capabilities of the MG5_aMC framework in order to have the possibility to compute NLO EW corrections and more in general QCD and EW complete-NLO predictions in the 4FS. This feature will be included in a future release of the code. To the best of our knowledge, this work is the first in which NLO EW corrections or complete-NLO predictions are computed in the 4FS for a complete process at hadron colliders. In the case of Hbb production, this work represents the first ever full computation of NLO EW corrections and complete-NLO predictions. While the impact of NLO EW corrections is found to be very small, once again complete-NLO predictions turns out to be much larger than naively expected values, due to the presence of new topologies, in this case the ZH and VBF ones.