NLO QCD Predictions for $Wb\bar b$ Production in Association with up to Three Light Jets at the LHC

In this article we present the next-to-leading order QCD predictions for $Wb\bar{b}+n$-jet ($n=0,1,2,3$) production at the Large Hadron Collider with $\sqrt{s}=13$ TeV. We work in the four-flavor number scheme with a non-vanishing bottom-quark mass and include all subprocesses at leading electroweak order as well as all heavy-fermion-loop effects. We show the impact of QCD corrections for total as well as differential cross sections and make an assessment of theoretical uncertainties of $Wb\bar{b}$ production viewed as an irreducible background to $H(\rightarrow b{\bar b})W$ studies. For the calculations we have employed an upgraded version of the BlackHat library which can handle massive fermions in combination with SHERPA. Our results can be explored through publicly available $n$-tuple sets.


I. INTRODUCTION
Particle phenomenology at the Large Hadron Collider (LHC) requires the measurement of proton-proton collisions with diverse final states including photons, leptons, heavy and light jets and missing transverse energy. In this article we provide theory predictions for W bb production in association with up to three light jets at the LHC in the Standard Model (SM) including next-to-leading order (NLO) QCD corrections. These processes have rich collider signatures, including almost all final-state objects mentioned above, and so they provide a natural test ground for precise measurements of complex signatures at the LHC.
Specifically, they appear as main reducible and irreducible backgrounds (together with nonresonant top-pair production) to HW associated production, with the Higgs boson decaying into a bottom-quark pair (bb). Currently this process receives increased attention by both LHC experiments, given its relevance for constraining the coupling of the Higgs to b quarks.
Providing results as a function of the light-jet multiplicity is particularly beneficial, as one expects that universal behavior appears at large multiplicities [1][2][3][4], as for example studied in vector-boson production in association with multiple light jets in ref. [5].
The first studies of NLO QCD predictions to W bb production appeared about twenty years ago [6] and were already included in the first version of the MCFM program [7], being important in particular for Higgs searches at the Tevatron. Those initial predictions were performed in the context of the four-flavor number scheme (4FNS) though they were performed in the approximation of massless quarks (employing the one-loop helicity amplitudes of ref. [8]). Results including full b-mass effects appeared first in refs. [9,10] with on-shell W bosons and subsequent refinements in refs. [11,12]. NLO QCD corrections for W bb + 1-jet production were obtained in ref. [13,14]. Also studies with more inclusive samples of b jets have been carried out. In ref. [15] NLO QCD corrections were computed for W + b + 1-light jet production, while NLO QCD results were presented for W production in association with a single b jet in refs. [16,17].
Experimental measurements have been carried out by the CDF experiment for W production with one or two b jets [18] (including samples with light jets), by the D0 experiment for inclusive W production with a single b jet [19], by the ATLAS experiment for W with up to two b jets [20] and also by the CMS experiment for W and two b jets [21,22]. For a review on these experimental measurements as well as theoretical predictions of vector boson production in association with b jets see ref. [23].
Here we present for the first time NLO QCD corrections to W bb+2-jet and W bb+3-jet production, and for completeness recompute the cases with zero and one light jet. These O(α s ) corrections are computed to the corresponding leading-order (LO) results at O(α 2+n s α 2 f ), with n the numbers of light jets in the process. At higher orders in the fine-structure constant α f , the same signatures can be obtained from processes involving top quarks which we do not consider here. Since the early calculations of inclusive W bb production [6,9,10] it has been observed that the NLO QCD corrections are quite large. This is mainly due to the opening of a gluon-initiated channel as part of the real contributions to the NLO QCD corrections, as well as to the release of a LO kinematical constraint which fixes the p T of the W boson to that of the bb system. Put differently, those processes suffer from giant K-factors [24]. The corrections to W bb + 1-jet production on the other hand show better behavior [13], and it is expected that for even larger light-jet multiplicities, the NLO QCD predictions will present more stable jet-scaling properties. We show that this is indeed the case in our results.
Attempts to obtain reliable predictions in spite of giant K-factors for W bb production have initially focused on exclusive analyses [9], including jet vetoes. But sensitivity to the p veto T cut tends to spoil the convergence of the perturbative series [25]. In this article we employ exclusive sums [26] as an alternative to stabilize these predictions, considering exclusive combinations of up to two light jets. The perturbative series for exclusive-sum observables is better behaved, as confirmed by comparison to LHC data for W +1-jet production [27,28].
Results based on exclusive sums are expected to actually contain some large next-to-nextto-leading order (NNLO) corrections, and so in the lack of a full NNLO QCD study of W bb production, our results represent a useful parton-level prediction for W bb observables. We put particular emphasis on observables associated to H(→ bb)W production, that is we study the p bb T , p W T , and M bb exclusive-sum distributions. We obtained our results using the BlackHat library [29] after a significant upgrade and the addition of new algorithms [30]. Most notably, the new version of the program allows to compute tree-level and one-loop matrix elements including massive fermions. This has been achieved by the implementation of the unitarity method [31][32][33] extended to massive particles [34], using a variant of the numerical unitarity approach [29,[35][36][37][38] that introduces a modified spectrum of particles in the loop. For the real-emission corrections we use the massive-dipole formalism [39], as implemented in the COMIX package [40] which is part of the SHERPA Monte Carlo program [41]. We store our results in a set of n-tuple files [42] which allow fast a-posteriori studies of our results, including scale variations, partondistribution functions (PDFs) and α s reweighting, and evaluation of additional infrared-safe observables.
This paper is organized as follows. In section II we give our calculational setup, details of the matrix-element computation as well as kinematical information, coupling schemes and other input parameters employed. In section III we present our results for total and differential cross sections and study their renormalization-and factorization-scale dependence.
Section IV presents a series of observables associated to HW production measurements and we show predictions of exclusive sums which combine cross sections of distinct light-jet multiplicities and assess their theoretical uncertainty based on scale dependence, higher-order contributions and PDF errors. We give our conclusions and outlook in section V.

II. CALCULATIONAL SETUP
We compute the W bb + n-jet (n = 0, 1, 2, 3) processes in NLO QCD, with the charged vector boson decayed to a massless charged lepton and its corresponding neutrino. The partonic subprocesses of Born and virtual contributions are obtained from the following subprocesses, for the seven-parton amplitudes are shown in fig. 1. The real-emission subprocesses are obtained from the above list by adding a gluon or replacing a gluon by a Q Q pair.
We obtain fixed-order parton-level predictions and do not include parton-shower effects.
All observables considered will be constructed from events containing exactly two observable b jets, defined in an infrared safe way [43]. We do not introduce any corrections due to possible mistagging of heavy and/or light jets.
Representative diagrams for two subprocesses contributing to pp → W bb+3-jet production.
The diagram (b) displays a contribution from closed loops of top and bottom quarks.
Below we give more technical details about our results including a brief description of the one-loop matrix-element computation in the newly setup BlackHat library, details on the renormalization schemes considered, numerical stability, validation, Monte-Carlo integration, input parameters, the considered observables as well our choices for the renormalization and factorization scales. We end this section with a brief assessment of b mass effects in our calculation.

A. Virtual Matrix Elements
A new version [30] of the BlackHat library [29] is used to compute the required virtual matrix elements, which includes significant upgrades for the computation of loop amplitudes with internal and external massive particles. The library uses the unitarity method [31][32][33] and its extension to massive particles [34] in order to compute loop amplitudes numerically.
These methods have already been applied by a number of groups for analytic as well as numerical computations of massive amplitudes (see e.g. [11,44,45]). The present implementation is based on the numerical unitarity approach [29,35,36,38] and its extension to massive quarks [37]. In addition, a prescription to map the higher-dimensional Dirac algebra into four-dimensional objects is used, which can be shown to be equivalent to the one given in ref. [46] with some modifications.
We apply the color decomposition of partial amplitudes in terms of primitive amplitudes [47,48]. Integrands A( ) of the primitive amplitudes are parameterized as [35,36,38], vanish upon integration over the loop momentum, and master terms, which are associated to master integrals. We customarily choose these terms such that the set of master integrals is formed by 1-, 2-, 3-, and 4-point scalar integrals, as well as by "extra" integrals originating from squares of (D − 4)-dimensional parts of the loop momentum [38,49].
Unitarity methods obtain the above numerators directly from on-shell tree amplitudes, thus avoiding an explicit reduction of tensor and higher-rank integrals, as in the associated OPP reduction method [35]. Suitable loop-momentum parameterizations allow to set n ≤ 5 loop propagators on-shell. On these on-shell phase spaces the amplitude factorizes into a product of n tree amplitudes such that in these limits the integrands in eq. (2.2) can be computed from simple tree-level input.
We have implemented loop-momentum parameterizations for on-shell conditions associated to arbitrary (real or complex) masses. We compute the required tree amplitudes numerically via Berends-Giele (BG) off-shell recurrence relations [50]. In particular, this allows for an efficient and flexible tree generation for both complex and D-dimensional momenta. Special attention needs to be paid to the computation of single cuts for the extraction of tadpole coefficients and double cuts for the extraction of coefficients of bubbles with a single on-shell leg. When computing the required single and double cuts, explicit divergences associated to self-energy insertions on external legs are encountered. These contributions are removed and accounted for by mass and wave-function renormalization. In order to do this, we follow ref. [37] and adjust the tree-diagram generation to remove these contributions from the cuts. Alternative approaches have been presented in ref. [51] and, recently, in ref. [52].
We compute intermediate expressions in the four-dimensional helicity (FDH) scheme.
The limit D s → 4 in eq. (2.2) has thus to be taken after the Dirac and vector algebra are evaluated. We address this with an approach based on D s -dimensional unitarity cuts [38].
Whenever possible we reduce the D s -dimensional algebra and states to lower dimensions. This is in close analogy to the known decomposition of the D s = 6 dimensional gluon amplitude to a D s = 5 amplitude plus a scalar contribution [32]. By doing so, we avoid an overhead of numerical computations in higher-dimensional representations. Our implementation can be equivalently derived from the four-dimensional (re-)formulation of FDH [46] with some modifications. Details will be discussed elsewhere [30].
Finally, we implemented the one-loop scalar integrals based on refs. [53][54][55]. The integrals are setup for evaluation in higher-precision arithmetics and have been systematically checked by comparing against the OneLOop Fortran library [56], as well as against the massless one-loop scalar integrals implemented in earlier versions of BlackHat.

B. Renormalization
We use the FDH variant of dimensional regularization in intermediate steps to regularize UV and IR divergences. At the end we convert the renormalized amplitude to the 't Hooft-Veltman (HV) scheme [57].
In table I we summarize the renormalization counterterms in the FDH scheme. For all external states we use the on-shell scheme and for the QCD coupling we use the M S scheme.
The gluon wave-function-renormalization counterterm receives contributions from all active heavy quarks. In the decoupling limit this counterterm is set to zero. In the 4FNS, the number of light flavors N f is set to four, and the number of heavy flavors N h is set to two, accounting for bottom-and top-quark loops. For each heavy quark included we also add a finite decoupling shift. Thus the full renormalized amplitude A (ren) is obtained from an amplitude with renormalized quark masses, A

Scheme Counterterm
Heavy quark wave function on-shell Light quark wave function on-shell 0 (UV+IR cancellation) Quark mass on-shell Gluon wave function on-shell  The renormalization of quark masses cannot be represented as a contribution proportional to the tree amplitude. We explicitly compute mass counterterm contributions using a dedicated recursive tree-like computation at the level of primitive loop amplitudes.
The conversion to the 't Hooft-Veltman scheme is performed by a finite shift [58] A (ren) 4) where N q stands for the number of light external quarks in the respective amplitude.

C. Numerical Stability
The numerical unitarity approach has proven to be numerically stable even when computing high-multiplicity one-loop matrix elements (see for example [59,60]  We produce histograms of the logarithmic relative error δ: for sample subprocesses of the most complex types in our calculations. The superscripts "prod" and "HP" mean normal-production evaluation and quadruple floating-point evaluation respectively. In figs. 2 and 3 we show a numerical stability plot for W bb + 2-jet and W bb + 3-jet production respectively. For the W bb + 2-jet process we include the two types of subprocesses: those associated to the four-quark (0 → W bbqq gg) and the six-quark (0 → W bbqq QQ) matrix elements, respectively. Similarly, for the W bb + 3-jet process we show results for the four-quark (0 → W bbqq ggg) and six-quark (0 → W bbqq QQg) matrix To control the precision of the computation we implemented a rescue system which identifies phase-space points which lead to numerically unreliable computations by performing checks at several stages of the calculation * . Whenever any of those checks fail we switch the computation to use higher-precision arithmetics locally, i.e. only for parts of the computation which failed the check. Typically, the fraction of time spent on these recomputations is small.
Overall we observe very good precision, with less than 1 in 10 4 phase-space points computed to an accuracy worse than three digits † . Some of the few points in the low-precision tail can be traced to the fact, that the scale set by the bottom mass m b and other scales of the problem s ij at sufficiently high energy are separated by several orders of magnitude. * With some minor adaptations we use the same implementation as described in [29,59]. † This precision target is chosen to allow for matrix element evaluation to a much better precision than the achievable statistical error in Monte Carlo integrations.
Ratios of the form (m 2 b /s ij ) k enter loop-momentum parameterizations and can cause loss of precision. We have observed that the contribution of the points in the low-precision tail is several orders of magnitude smaller than the total statistical errors for the observables studied.

D. Validation
We have carried out a number of checks with the new implementation of BlackHat [30].
In particular, we have systematically reproduced all massless results that were carried out with all the earlier versions of the library, and found excellent agreement.
We also apply a number of internal checks. For each primitive amplitude we check consistency of the integrand reduction by attempting to compute a high-rank tensor coefficient, which is known to vanish from power counting. We also check that the matrix elements reproduce the known IR and UV singularity structure [62]. Both of these checks are also used to control the precision of the computation. Explicit cancellation of the infrared poles of the renormalized one-loop matrix elements is also checked by comparing to the integrated subtraction terms, after PDF renormalization, as performed with the SHERPA library.
We also reproduced the results presented by Ellis et al. in [37] at the level of primitive amplitudes. Even more, we have cross checked fully interfered, matrix-element squared results against publicly available matrix-element generators. This includes comparisons of all one-loop matrix elements necessary for pp → tt + (≤ 2)−jet and pp → bb + (≤ 2)−jet at NLO QCD against the Recola [63] and OpenLoops [64] libraries (both using the Collier library [65]), as well as all the ones needed for W bb + (≤ 3)−jet with Recola.
Finally, we made dedicated comparisons at the level of physical observables against the MCFM program [7] for the inclusive production of W bb at NLO QCD at the LHC with √ s = 13 TeV. Agreement was found at the permil level for total cross sections and differential distributions.

E. Monte Carlo Integration
NLO QCD corrections require the integration of Born, virtual and real-radiation matrix elements. We provide virtual matrix elements to the SHERPA Monte Carlo program [41], which are then combined with Born and real-radiation matrix elements and integrated consistently.
The results in this article were obtained by exploiting the color structure of the matrix elements to increase the efficiency of the phase-space integration. To this end we split up the virtual matrix elements into leading and subleading color terms [48,59] and sample them independently over phase space. Given the color suppression of the subleading color terms these computationally expensive contributions need to be evaluated less often to reach a given total statistical-integration error.
The infrared singularities of the real-emission corrections have to be cancelled explicitly against their counterparts in the virtual matrix elements before numerical integration. We use the subtraction scheme based on massive dipoles [39] in our calculation. In particular, we employ the automated implementation in the COMIX package [40], which is part of the SHERPA library [41]. The latter implementation has been checked extensively in the recent calculation of tt production in association with up to three light jets [66].
We also employ the SHERPA library to handle the subprocess generation and mapping.
Furthermore, we rely on SHERPA for the phase-space integration and use its internal Analysis package. We compute fixed-order parton-level predictions and include neither partonshower effects nor hadronization corrections or other non-perturbative effects. We store our results in flexible ROOT [67] n-tuple files, which allow for a-posteriori variations of the strong coupling, PDF sets and choices of renormalization and factorization scales of the NLO QCD results. The format was developed in ref. [42].
We evolve α s (µ) with the QCD beta function for four massless quark flavors, for all µ. This is achieved by introducing a decoupling shift for massive quarks (see above). We use a one-loop running of α s at LO and a two-loop running at NLO. We choose m b = 4.75 GeV, consistent with the PDF sets, and set the top-quark mass in closed-loop contributions to m t = 172 GeV.  [71].

Parameter
Value We work at leading order in the electroweak coupling and fix the W -boson couplings to fermions with the SM input parameters as specified in table II. We use the G µ scheme [70] and compute the parameters α f (M Z ), sin 2 (θ W ) and g 2 W using the tree-level relations The lepton-pair invariant mass follows a relativistic Breit-Wigner distribution, with the mass and width as specified in table II. We approximate the Cabibbo-Kobayash-Maskawa (CKM) matrix by a unit matrix. This results in a small change of the total cross sections for the setup we use, as estimated by LO evaluations with the full CKM matrix. Indeed we find that these differences are of the order 1% for W bb production and below 0.5% for W bb + 1-jet and W bb + 2-jet production.
All light quarks (u, d, s, c) are treated as massless. We do include contributions from virtual bottom and top quarks, and we confirm the expected percent-level effect on crosssections [72][73][74]. We work with a single massless lepton pair, an approximation that can be applied to the electron or muon families.

G. Kinematics, Observables and Exclusive Sums
For completeness, we state the definitions of the standard observables used in our analysis.
The pseudorapidity η is given by with the polar angle θ with respect to the beam axis. The angular separation between any two objects (partons, jets or leptons) is denoted by where ∆φ is the difference in the azimuthal angle in the plane transverse to the beam axis and ∆η the difference in the pseudorapidities. The total partonic transverse energy is given where the sum runs over all final state partons j, independent of whether or not they are inside a jet that passes the cuts. The scalar transverse momentum p T of a parton is given by p T = p 2 x + p 2 y and the transverse energy of the W boson is computed from (2.10) The total partonic transverse energy is not directly measurable. Nevertheless, it is a suitable candidate for choosing unphysical renormalization and factorization scales since changing the cuts does not affect the matrix element at a given phase-space point. Jet invariant masses are defined by where jets are labeled in order of decreasing transverse momentum p T . The transverse mass of the W boson is computed from the kinematics of its decay products In section IV, we show predictions based on exclusive sums [26], which combine exclusive (σ exc n ) and an inclusive prediction (σ inc n ) for distinct light-jet multiplicities n. We study exclusive sums for observables in W bb production (0 light jets) and use our high-multiplicity results to define the exclusive sums labeled 'NLO+' and 'NLO++', (2.13)

H. Dynamical Scale Choice
We present renormalization and factorization scale dependence of cross sections using correlated variations of a central scale choice µ 0 by factors of (1/2, 1/ √ 2, 1, √ 2, 2). We have also explored independent variations of those scales and find similar results as with the correlated variations. We chooseĤ T /2 (2.9) as the functional form for the central scale µ 0 and keep factorization and renormalization scales equal, µ R = µ F = µ 0 . The scale µ 0 =Ĥ T /2 has proven to be a sensible choice as it tends to reduce shape changes and the global size of quantum corrections when going from leading to next-to-leading order (see for example [5,60,73]). In general, NLO corrections are less sensitive to the choices of scale, as long as the scale reflects the hardness of the Born process [59]. We will label results for leading-order QCD with the central scale µ 0 =Ĥ T /2 by 'LO' and the corresponding results at next-to-leading order QCD by 'NLO'.

I. Effects of a Finite b Mass
Mass effects in W bb have been studied since the early NLO QCD calculations in ref. [9].
They are expected to be small when two well-defined b jets are considered, and ratios of m 2 b to typical invariants are small. Nevertheless, their contributions are fundamental when studying inclusive b-jet production at hadron colliders (see for example refs. [16,17]).
In order to highlight these effects, in fig. 4 we show the ratio of a computation performed in the 4FNS consistently keeping the mass of the b quarks, to that of a corresponding massless calculation performed with massless b quarks in the five-flavor number scheme (5FNS). In the latter we use the PDFs from CT14 [68], denoted by CT14llo at LO and CT14nlo at NLO. We notice that for values of M bb above 50 GeV, the ratios stabilize rapidly at about 0.95 for W bb production and at 0.9 for W bb + 1-jet production, while for values below we have a strong deviation with the massless calculation more than doubling the massive one. This is to be expected as phase space constrains the production of massive b quarks in these regions and also m 2 b /M 2 bb terms in the matrix elements can be important. The mass effects are stable with respect to quantum corrections, as we can deduce from the similarity of the LO and NLO ratios. We notice that the deviation from 1 at large M bb is smaller than the scale-dependence bands of the NLO results, which can be used as a proxy of unaccounted higher-order corrections.
The observed behavior is very similar at what was studied in ref. [9] in the case of W bb production, and here we extend it to W bb + 1-jet production. It is important to mention that although the computations in fig. 4 are consistent results in the 4FNS and in the 5FNS, they have the same diagrammatic content at Born level, and in particular there is no binitiated subprocess in the massless calculation. For that reason the comparison presented is attributable to b-mass effects in the matrix elements and in phase-space generation, together with the corresponding differences from the PDFs and their corresponding running couplings.
A more systematic 4FNS vs. 5FNS comparison including W bb + 2-jet and W bb + 3-jet production would be more relevant to compare the two schemes, which we leave to future work. A future study of our results for more inclusive samples of W production in association with b jets will also be interesting (in this study we focus in signatures with exactly 2 b jets).

III. RESULTS FOR W bb PRODUCTION IN ASSOCIATION WITH LIGHT JETS
In this section we present NLO QCD results for W bb + n-jet (n = 0, 1, 2, 3) production in pp collisions at √ s = 13 TeV, the experimental configuration of the LHC Run-II. We present results for a set of distributions and apply the following cuts The cuts are applied to both light and b jets. The renormalization and factorization scales are chosen to be equal and set on an event-by-event basis by µ =Ĥ T /2, according to eq. (2.9). We define our jets by employing the anti-k T jet algorithm [75] with R = 0.4, as implemented in the FastJet package [76].

A. Total cross sections
In table III, we present total partonic cross sections, employing the kinematical cuts of eq. (3.1), for inclusive production of both W − and W + accompanied by two b jets and zero to three light jets. The numerical integration uncertainty is given in parenthesis and the scale dependence is quoted in superscripts and subscripts. We also show the ratio of NLO over LO results, so called K-factors, in separate columns.
LO cross sections display a large scale sensitivity, reaching up to 60% for W bb + 3-jet production. We note that the scale dependence of the LO cross section for W bb is around 20% while the NLO QCD corrections increase the total cross section by a factor of 2. This clearly highlights that scale dependence is in general not representative of the associated theoretical uncertainties. In this case, the large quantum corrections can be understood as a result of the opening of gluon-initiated channels [6,9,10]. Also for W bb + 1-jet a gluongluon initiated channel is opened up, but with milder impact, and for the larger multiplicity processes all subprocesses are present at LO. Hence, quantum corrections are milder for these processes. Furthermore, kinematical constraints at LO are only present for W bb production, as we will discuss for example for the p bb T and p W T observables in section IV. As a consequence, we expect quantum corrections for processes with even more light jets to be under relatively good pertubative control.  The last four columns give jet production ratios for both W − bb as well as W + bb in association with n light jets. These ratios are taken for the cross section of a given process to that with one less jet. The number in parenthesis gives the corresponding statistical integration error.
In table IV we  Finally we also explore in table IV the jet ratios in W ± bb production in association with light jets. Similarly to studies of these ratios in W +n-jet (light jet) production [5], we observe that the results for n = 1 are special given the large NLO corrections for W bb production.
The opening of an initial-state channel makes the W bb + 1-jet/W bb ratio clearly sensible to quantum corrections. In the light jet study [5] this was the case for the ratio (W + 2jet)/(W + 1−jet), and a full study of jet-ratio universality needed the completion of the NLO QCD correction to W + 5-jet production [60]. Similarly, in W bb inclusive production, it might be interesting to explore the NLO QCD corrections to W + bb + 4-jet production in the future.

B. Scale Dependence
In fig. 5 we study the dependence of total cross sections in W − bb and W + bb production in association with up to 3 light jets on the renormalization and factorization scale. We employ the central dynamical scale µ 0 = µ R = µ F =Ĥ T /2. The scale variations observed for W + and for W − are very similar. The LO cross sections have a monotonically increasing scale dependence, for n ≥ 1. As we observed in the previous subsection, the scale dependence of W bb production is special.
We choose the dynamical scaleĤ T /2 which on average increases monotonically with multiplicity. For vector boson production in association with massless jets this scale choice has been observed to produce stable NLO results over a wide range of kinematical configurations relevant to the LHC and future colliders [60,72,78,79]. For the LHC in particular, it has been observed that for massless jet production the scaleĤ T /2 typically produced NLO cross sections lying on the locus of the scale-dependence curves. Here we observe that for W bb + n-jet (n = 0, 1, 2, 3) production, the NLO cross section at the central scale appears consistently on the right of the scale-dependence plateau. We can assert, in particular considering the similarities of the massive and massless results studied in section II I, that this difference has little to do with the presence of a massive jet, and it is actually due to the dominant type of subprocess. For light-jet production those are the ones with a single quark line, while in the case of W bb + n-jet (n = 0, 1, 2, 3) production the dominant subprocess are those with two quark lines (those are the subprocesses with most gluons allowed).
Another interesting difference between W production in association with light jets and W bb production with multiple light jets, is that for the former the leading-color approximation for one-loop matrix elements gave a very good approximation for physical observables (at the level of 1 to 3%). Contrary to that, W bb production with light jets is largely dominated by virtual contributions in our setup, and so the leading-color approximation is at the order of 10% for physical observables. That is why all of our results in this article include full-color information, and we only exploit the decomposition in a color expansion for efficiency of the computation. We again attribute this difference to the unlike dominant subprocesses. In the bottom panels we show the scale-dependence bands normalized to the NLO result, in blue for LO and dark gray for NLO.

C. Differential distributions
In this section, we describe NLO results for several differential distributions and thereby analyze the impact that quantum corrections have on fixed-order predictions over phase space. We generally show results only for one of the W ± charges, as the structure of the corrections are similar between them.
In figs. 6 and 7 we show the jet-p T spectra of the leading and subleading b jets (ordered by p T ) respectively, for inclusive W − bb production in association with n = 0, 1, 2, 3 jets. The NLO corrections show quite some structure beyond the K-factors studied at the level of the total cross sections in the previous subsection. We observe shape differences in most of the p T distributions of the b jets, in a way that make the LO predictions usually harder (with the exception of the leading b jet p T in W bb production). Nevertheless, the LO/NLO shape difference is clearly reduced for the process with highest multiplicity, W bb + 3-jet production, a feature that shows up persistently in the following observables. We notice that the scale dependence of the NLO results is reduced compared to the LO results (apart from W bb, as discussed for total cross sections). In the high multiplicity samples, the NLO results lie inside the LO bands.
In fig. 8 we show the p T distributions of the softest light jet in inclusive W + bb+1, 2, 3-jet An interesting observable for many scenarios of physics beyond the SM (BSM), as well as for experimental studies at hadron colliders, is that of the total hadronic activity in a detector. In fig. 9 we show the distribution in this observable, including all hadronic activity from the light and b jets in W − bb+0, 1, 2, 3-jet production. Large and phase-space dependent NLO corrections appear for W bb as we would expect from previous discussions. Interestingly, in W bb + 1-jet production a remnant of these large effects appears in this observable. The corrections are not as large as for the former, but still at around 1 TeV for H jets T , we see a differential K-factor reaching two, though the shape difference seems to end at about 400 GeV. The large-multiplicity processes on the other hand show much less structure, related to the kinematically unconstrained nature of their LO configurations, which contain any W soft enhancements starting at LO.
Finally, to end this section, we show in fig. 10 the distribution on the ∆R separation between the first b jet and charged lepton l − . Most of the angular variables that we have studied are similar to this one, which shows little structure in the QCD corrections. We only find effects when a certain kinematic constraint is imposed at LO and released by the corrections, as it is the case on the left most plots of fig. 10. In the case of the ∆R bl − at LO in W − bb production, the parent W and gluon that give rise to the leptons and b jets are produced with ∆φ (the difference in azimuthal angle) equal to π. Also, the ∆η distribution peaks at around zero and decreases monotonically. The resulting ∆R bl − distribution thus has the feature of a sharply decaying distribution at LO. All those constrains are lifted by the real corrections and do not appear at all in W bb+1, 2, 3-jet production.

IV. BACKGROUNDS TO HW PRODUCTION
So far all LHC measurements of Higgs-boson properties appear in good agreement with predictions of the SM (see for example ref. [80]). One of the properties that will be important to constrain further is the coupling strength of the Higgs boson to b quarks. Given that a Higgs boson with mass M h around 125 GeV is supposed to decay more than half of the time into a bb pair, it is of great importance to constrain the Yukawa coupling y b and consequently learn about the Higgs boson's total decay width. In the main production channel of the Higgs, through gluon-gluon fusion, one faces a large background from pure QCD bb production. Therefore, considering the associated W bb production gives an extra handle to detect the Higgs decaying to b quarks (see the recent measurement by ATLAS in ref. [81]). This of course, as long as the irreducible backgrounds for W bb production in the SM can be kept under control. The predictions provided in this article aim to contribute to these studies.
Some of the key observables for HW analyses are those associated to the bb system, in particular p bb T and M bb . When producing a Higgs, those are associated with the p T distribution and the invariant mass of the Higgs. In addition, distributions that help to characterize the accompanying W boson are important, for example p W T . In this section we study those three observables with our high-multiplicity results. At high energies the presented spectra are enhanced with resonant top contributions (see for example the recent study [82]). Nevertheless, in the context of HW production the focus is on the non-resonant backgrounds. The non-resonant top contributions are sizable, and can be of similar order to the ones presented here. We leave studies of non-resonant top contributions to future work.
The NLO QCD correction to W bb production have large contributions associated to processes with an extra light jet [6,9,10]. A way to handle those contributions was to obtain exclusive results in the number of light jets [9], which explicitly vetoed events with extra jets, but that prescription suffers from its sensitivity to the p veto T cut [25]. In this article we consider a different approach, using the 'exclusive sums' technique [26]. Instead of imposing a veto cut to stabilize the predictions, this approach replaces the extra lightjet contributions to generic observables, which are effectively LO, by their corresponding results including NLO QCD corrections. In higher-order corrections such contributions are naturally added. However, as they are hard to obtain, we use the above approximation and analyze the impact in our predictions.
The 'exclusive sums' technique is expected to give improved predictions when tree-like contribution, with an extra light jet, to NLO corrections are large. Notice that in measurements of W +light jets some of the predictions from exclusive sums have been compared to LHC data, see for example [27,28], usually in the context of W + 1-jet production. By now those computations are outdated, given the recent NNLO QCD calculation presented in ref. [83]. It is important to mention that for the comparison to LHC data the application of a parton-shower can be studied [13], or using a matched and merged version for example with the MEPS@NLO [84] or FxFx technique [85].
We will focus on predictions for W bb+X production. Fixed-order results for those will be denoted as usual with the labels LO and NLO. The exclusive sums we employ are defined in dσ/dp bb T [pb/GeV] Eq. (2.13), for which we will use the labels NLO+ and NLO++. We will use the latter only as a proxy for the size of even higher-order corrections, that is as an estimator of theoretical uncertainty. Our main predictions will be that of NLO+.
For completeness, to characterize theoretical uncertainties, we will also explore the PDF uncertainty associated to the observables under consideration, even though they turn out to be subleading. In order to estimate the PDF uncertainty we use the error sets from the pseudo-PDF set PDF4LHC15 nlo nf4 30 [86]. Given the smallness of the PDF errors compared to other theoretical uncertainties, we do not go beyond this restricted error set for our estimations. Other sources of uncertainties are the values of m b and α s , but those are expected to be rather small, e.g. when compared to missing higher-order corrections.   In fig. 11 show a difference in shape, particularly at low values of p bb T . We find that the higher-order corrections estimated through scale variations and by the NLO++ results are of the same order of magnitude, at the level of 10%.
For the p bb T observable, one obseves the release of kinematical constraints at NLO (which tights up p bb T and p W T ). Since at LO the massive b-quark pair originates in a gluon splitting, the kinematical constraints in W bb are similar to those appearing in V + 1 light jet (see e.g. refs. [60,72,78]). Real-radiation emission relaxes this constraint and yields large corrections at NLO through a soft enhancement, which gives rise to a giant K-factor [24]. These characteristics of the NLO results for W bb production are then expected [87] and require resummation although fixed order corrections are known to improve the predictions [88].
In fig. 12, we show in the same manner the distribution in the transverse momentum of the vector boson p W T . The results are again similar to what we encounter for p bb T , with the NLO+ uncertainty estimation marginally overlapping with the NLO predictions and with its scale sensitivity of the order of the NLO++ predictions. The estimation of theoretical uncertainties is about 17% in the range of p W T shown (as compared to the 25% of the NLO results). Again, the PDF uncertainties are subleading, appearing at 3% and below.
Finally in fig. 13 we show the distributions in the M bb observable, which exhibit similar features to the observables studied previously in this section. In this case, the uncertainties associated with scale sensitivity and missing higher-order effects appear at about 10% or slightly smaller. The PDF uncertainties are tiny, being at 1% or smaller.

V. CONCLUSIONS
In this paper we presented the first NLO QCD results for W bb + 2-jet and W bb + 3jet production. For completeness we have also included results with zero or one extra light jet at NLO QCD. These are the first results obtained with a new version [30] of the BlackHat library [29] which has been constructed to handle massive particles in one-loop   matrix elements.
We have found that NLO QCD corrections for W bb+2-jet and W bb+3-jet production are mild, with the associated K-factors dropping from 1.4 for W bb + 1-jet, to 1.2 for W bb + 2-jet, and 1.15 for W bb + 3-jet production. For the latter process the remaining renormalization and factorization scale dependence amounts to 60% at LO and 20% at NLO for the typical variation in a correlated or uncorrelated way of the scales by a factor of two. In contrast to this the K-factor for NLO QCD corrections is large in inclusive W bb production [6,9,10].
This difference led us to explore possible improvements to observables for the W bb process by using exclusive sums, in particular in the context of H(→ bb)W associated production.
We showed that exclusive-sum predictions for key observables give uncertainties associated to missing higher-order corrections at the level of 10% to 17%. We found that uncertainties associated to PDFs are subleading, at the level of 2% in the generic kinematical regimes studied.
The ATLAS and CMS collaborations have performed systematic measurements of vectorboson production in association with light jets (see for example refs. [27,89,90]) and these measurements have shed light on the precision of theoretical and experimental tools to describe processes with many objects in the final state. In the future, it will be interesting to use high-multiplicity processes including b jets, like the ones studied in this article, to verify both theory predictions and experimental techniques. We provided n-tuple sets for the predictions obtained by BlackHat and SHERPA for guidance in future analyses of the W bb + n-jet signatures.