$t\bar{t}b\bar{b}$ at the LHC: On the size of corrections and $b$-jet definitions

We report on the calculation of the next-to-leading order QCD corrections to the production of a $t\bar{t}$ pair in association with two heavy-flavour jets. We concentrate on the di-lepton $t\bar{t}$ decay channel at the LHC with $\sqrt{s}=13$ TeV. The computation is based on $pp \to e^+ \nu_e\, \mu^-\bar{\nu}_\mu\, b\bar{b} \,b\bar{b}$ matrix elements and includes all resonant and non-resonant diagrams, interferences and off-shell effects of the top quark and the $W$ gauge boson. As it is customary for such studies, results are presented in the form of inclusive and differential fiducial cross sections. We extensively investigate the dependence of our results upon variation of renormalisation and factorisation scales and parton distribution functions in the quest for an accurate estimate of the theoretical uncertainties. We additionally study the impact of the contributions induced by the bottom-quark parton density. Results presented here are particularly relevant for measurements of $t\bar{t}H(H\to b\bar{b})$ and the determination of the Higgs coupling to the top quark. In addition, they might be used for precise measurements of the top-quark fiducial cross sections and to investigate top-quark decay modelling at the LHC.


Introduction
The discovery of the Higgs boson at the LHC was only the start of the wide program, the main purpose of which is to identify the properties and couplings of this recently discovered particle. In the Standard Model (SM) the Higgs boson couples to the fundamental fermions via the Yukawa interaction with a coupling strength that is proportional to the fermion mass. Probing the coupling of the Higgs boson to the top quark, the heaviest observed particle in the SM, comprises a crucial test of the consistency of the Higgs sector. Furthermore, the Top-Yukawa coupling, denoted as Y t , might be used to constrain various models of physics beyond the SM (BSM) that very often predict a different coupling strength than the SM one. The latter is expected to be close to unity. Indirect constraints on the coupling between the top quark and the Higgs boson are available from processes including virtual top quark loops. Here the best example comprises Higgs boson production through gg fusion. On the other hand, the Y t coupling can be probed directly in the associated production of the Higgs boson with the tt pair, the process which has recently been observed by both the ATLAS and CMS collaborations [1,2]. Even though pp → ttH contributes only about 1% to the total pp → H production cross section, it offers a very distinctive signature. For the Higgs boson with the observed mass value the dominant decay mode is H → bb with the branching ratio of 58% [3]. This decay mode is additionally sensitive to the coupling of the Higgs boson to the bottom quark but it is not easily accessible experimentally. Nevertheless, both ATLAS and CMS reported searches for ttH production in the bb decay channel of the Higgs boson [4,5]. The main experimental challenge for this channel is the correct identification of the candidates for the Higgs boson decay from the so-called combinatorial background. The latter is responsible for a substantial smearing of the Higgs boson peak in the bb invariant mass spectrum. Further challenges include the possibility of misidentification of light jets with b-jets and problems with the control of various SM backgrounds, see Figure 1 for examples of Feynman diagrams for ttH(H → bb) production and processes that lead to the same ttbb final state 1 . With the help of b-jet tagging algorithms as well as boosted top quarks and Higgs boson it is possible to isolate the contribution of the ttH(H → bb) process from the most general reducible background represented by ttjj production and from the irreducible Z-peak background [7]. Nevertheless, the pp → ttH → ttbb process suffers greatly from the ttbb background, that is the most important irreducible background process for this SM Higgs boson channel. In addition, searches for four top-quark production (tttt) are also affected by the QCD ttbb background [8,9]. Consequently, measurements of ttH(H → bb) and tttt at the LHC would benefit considerably from a better understanding of the QCD ttbb production process and from the improved modelling of top-quark decays.
The pp → ttbb process is also interesting from the theoretical point of view. This is due to the presence of two very different and distinctive scales, the top-quark mass m t and the b-jet transverse momentum. The former governs tt production and subsequent top-quark decays, the latter describes the two b-jets coming from the g → bb splitting. However, this rather straightforward picture is not adequate anymore once the contributions from off-shell top quarks and W gauge bosons are consistently taken into account. Away from the tt threshold and for large values of M (bb) various other mechanisms start to play a non negligible part, that makes the ttbb process truly multi-scale production.
Calculations at NLO in QCD are available for stable top quarks already for some time [10][11][12][13][14][15]. They suffer from large uncertainties in the choice of factorisation and renormalisation scales that are of the order of 33% and large NLO corrections of the order of 77%. The latter is mainly due to the gg initial state. Indeed, NLO QCD corrections to the subprocess initiated by qq annihilation only are small with the K-factor of the order of K = σ NLO /σ LO ≈ 2.5% whereas the NLO theoretical error due to scale variations is at the 17% level [10]. Unfortunately, the large scale variation and the size of the corrections themselves for the full pp process imply that a full NNLO study would be indispensable. As the latter will remain out of reach in the nearest future, it seems that, as already suggested in Ref. [13], additional kinematic restrictions to the bb system, e.g. M (bb) 100 GeV or p T (bb) 200 GeV as motivated by the studies of the SM Higgs boson with the mass of m H = 125 GeV, must be introduced in order to reduce large theoretical uncertainties and higher-order QCD corrections. Alternatively, a veto on extra jet radiation might be carried out to achieve the same goals. In Ref. [16] even calculations of NLO QCD corrections to the pp → ttbbj process were presented. For integrated NLO cross sections scale uncertainties at the level of 25% were obtained, whereas, the usual K-factor was calculated to be 1. 45. Having an additional resolved jet present allowed to investigate in detail the modelling of recoil effects in ttbb production. We note here that, generally speaking, the pp → ttbbj process at NLO entails information that can be used to gain some insights into the perturbative convergence of the inclusive ttbb cross section beyond NLO.
Computing higher-order QCD corrections to processes with stable top quarks, no matter how technically complex they might be, can only give us a general idea of the size of the NLO corrections. Such theoretical predictions cannot provide a reliable description of the top-quark decays and are not sufficient to detail QCD radiation pattern for this process. Thus, for more realistic studies radiative top-quark decays are needed. The first step in this direction was achieved, for the ttbb production process, by matching the ttbb matrix elements, calculated with massless or massive b quarks at NLO in QCD, to parton-shower (PS) programs [17][18][19][20]. In these studies, however, top-quark decays have been either completely omitted or performed within the Pythia MC framework. Either way, these predictions did not include tt spin correlations. On the other hand, in Ref. [21] NLO matrix element calculations matched to a parton-shower (NLO+PS) simulations of pp → ttbb were presented in the four-flavour scheme for stable top quarks and with (LO) spin correlated top-quark decays. In general, it was shown that matching and shower uncertainties were very small for this process once m b = 0 was considered. Moreover, neither K-factor nor the size of scale uncertainties was substantially changed by PS effects. Actually, the scale dependence uncertainties even increased slightly, reaching 40%. The experimental measurements of pp production cross sections for ttbb have been carried out by both ATLAS and CMS [22][23][24]. The measured inclusive fiducial cross sections generally exceed theoretical predictions for ttbb as provided by the already mentioned NLO + PS simulations. Nevertheless, they are compatible within the total uncertainties. Very recently, the computation of NLO QCD corrections to ttbb production in the di-lepton top-quark decay channel was reported in Ref. [25]. Specifically, higher order α s corrections to the e + ν e µ −ν µ bb bb final state at O(α 4 s α 4 ) were calculated for the LHC energy of √ s = 13 TeV. In Ref. [25] they have shown that at the level of cross section the NLO QCD corrections to ttbb were close to 100%. Furthermore, at the differential level on top of this overall shift, moderate shape distortions up to 25% were obtained.
The goal of the paper is manifold. Due to the complexity of the process we will present an independent computation of the complete NLO QCD corrections to the off-shell production of ttbb in the di-lepton top-quark decay channel. In the first step, we use the same scale choice and the SM input parameters as well as phase-space cuts to confirm the results presented in Ref. [25] at the integrated and at the differential level. Our second goal is to extend this analysis to other selected renormalisation and factorisation scale choices and to different PDF sets. By using the error PDF sets we will additionally study the internal NLO PDF uncertainties. Furthermore, a stability test of LO and NLO fiducial cross sections with respect to the b-jet transverse momentum cut will be performed. Afterwards additional cuts will be introduced to examine their impact on the size of the K-factor and the theoretical uncertainties due to scale dependence. Another goal of the paper is to investigate the initial state bottom-quark contributions and their impact on the integrated and differential ttbb cross sections. To this end two approaches, the so called charge-blind and charge-aware heavy-flavour jet tagging, will be introduced.
The paper is organised as follows. In Section 2 we briefly describe the Helac-Nlo framework that we use for the calculation and discuss the cross-checks that have been performed. Our theoretical setup for LO and NLO QCD results is given in Section 3. Phenomenological results for the integrated and differential fiducial ttbb cross sections are discussed in detail in Section 4 and Section 5. The initial state bottom-quark contribution is examined in Section 6. The comparison with theoretical predictions presented in Ref. [25] is performed in Section 7. Finally, in Section 8 our results for the ttbb production process are summarised.

Description of the calculation and its validation
We consider the fully realistic e + ν e µ −ν µ bb bb + X final state. We consistently take into account resonant and non-resonant top-quark contributions and all interference effects among them. In addition, non-resonant and off-shell effects due to the finite W gauge boson width are included. Due to their insignificance we neglect flavor mixing. A few examples of Feynman diagrams contributing to the leading order process at O(α 4 s α 4 ) are presented in Figure 2. They are shown for the dominant gg partonic subprocess. NLO QCD corrections are calculated with the help of the Helac-Nlo Monte Carlo (MC) program [26]. This is the first computation of a 2 → 6 process (the decay products of the W 's are not counted, because they do not couple to colour charged states) carried out within this framework. Even though Helac-Nlo has already been employed for the calculations of NLO QCD corrections to tt+X, X = j, γ, Z(→ νν), W ± (→ ν) with full top quark off-shell effects included [27][28][29][30][31], these processes were at most 2 → 5 processes from the QCD point of view. For the gg → e + ν e µ −ν µ bb bb partonic reaction there are 3904 LO diagrams. For each qq → e + ν e µ −ν µ bb bb subprocess, where q stands for u d, c, s, we have 930 LO diagrams. The calculation of the LO scattering amplitudes is performed within the Helac-Dipoles package [32]. The phase-space integration is performed and optimised with the help of Parni [33] and Kaleu [34]. The produced top quarks are unstable particles, thus, the inclusion of their decays is performed in the complex-mass scheme [35][36][37][38]. It fully respects gauge invariance and is straightforward to apply. The resonant electroweak vector bosons are also treated in the complex-mass scheme.
We compute the virtual corrections using Helac-1Loop [39] and CutTools [40]. Specifically, one-loop amplitudes are generated with Helac-1Loop and are further reduced at the integrand level using the so-called OPP method [41] as implemented in CutTools. The most complicated one-loop diagrams in our calculations are octagons (8-point integrals). In the gg channel they involve tensor integrals up to rank six. In Table 1 the number of one-loop Feynman diagrams, that corresponds to each type of correction for the dominant gg partonic subprocess, is provided (see examples in Figure 3). They are obtained with the help of the Qgraf program [42] as Helac-Nlo does not employ Feynman diagrams for the amplitude calculations. We have cross-checked our results with the publicly available general-purpose MC program MadGraph5 − aMC@Nlo [43]. Specifically, we compared numerical values of the one-loop virtual corrections for a few phase-space points for gg and qq (q = u, d, b) partonic subprocesses. Virtual corrections come from the interference of the tree-level and one-loop amplitudes, summed over all colors and spins, and for N F = 5 massless quark flavours. We compared the finite parts along with the coefficients of the poles in . Additionally, coefficients for color and spin summed results for the I-operator were cross checked between Helac-1Loop and Helac-Dipoles. We have found perfect agreement in each case. At the one loop level the appearance of a non-zero topquark width in the propagator requires the evaluation of scalar integrals with complex masses, which is supported by the OneLOop program [44], used for the evaluation of the integrals. The preservation of gauge symmetries by this approach is explicitly checked by studying the Ward identities up to the one-loop level. For gg subprocess we perform this test for every phase-space point. Quadruple precision is used to recompute events which fail the gauge invariance check. For the qq subprocess we use the so-called scale test [45], which is based on momentum rescaling. Since we know how the recalculated amplitude scales when the momenta are rescaled, it is possible to compare the two results. Also in this case the test is performed for each phase-space point. For the failed points the amplitude is recomputed using higher precision. In addition, reweighting techniques, helicity and colour sampling methods are used in order to optimise the performance of the Helac-Nlo framework.
To compute the real corrections we isolate the singularities from the soft or the collinear parton emissions via subtraction methods for NLO QCD calculations. We employ Helac-Dipoles, which implements the dipole formalism of Catani and Seymour [46,47] for arbitrary helicity eigenstates and colour configurations of the external partons and the Nagy-Soper subtraction scheme [48], which makes use of random polarisation and colour sampling of the external partons. Two independent One-loop correction type Number of Feynman diagrams Self-energy 93452 Vertex 88164 Box-type 49000 Pentagon-type 25876 Hexagon-type 11372 Heptagon-type 3328 Octagon-type 336 Total number 271528 qq → e + ν e µ −ν µ bb bb g 9576 50 10 gq → e + ν e µ −ν µ bb bb q 9576 50 10 gq → e + ν e µ −ν µ bb bbq 9576 50 10 Table 2. The list of partonic subprocesses contributing to the subtracted real emission at O(α 5 s α 4 ) for the pp → e + ν e µ −ν µ bb bb + X process where q = u, d, c, s. Also shown are the number of Feynman diagrams, as well as the number of Catani-Seymour and Nagy-Soper subtraction terms that correspond to these partonic subprocesses.
subtraction schemes allow us to cross check the correctness of the real emission part of the calculation in an even more robust way. We use a phase-space restriction on the contribution of the subtraction terms as originally proposed in Ref. [12,49,50] for the Catani-Seymour scheme and in Ref. [51] for the Nagy-Soper one. We consider two extreme choices to cross check the independence of the results on this parameter. All partonic subprocesses that are taken into account for the real emission contributions are listed in Table 2, together with the number of the corresponding Feynman diagrams, the number of Catani-Seymour dipoles and Nagy-Soper subtraction terms. In each case, there are five times fewer terms in the Nagy-Soper subtraction scheme compared to the Catani-Seymour approach. The difference corresponds to the total number of possible spectators in the process under scrutiny, which are relevant only in the Catani-Seymour case.
Our theoretical predictions are stored in the form of modified Les Houches Event Files [52] and ROOT Ntuples [53]. Inspired by the ideas proposed in Ref. [54] each "event" is stored with sup-plementary matrix element and PDF information. Ntuples contain unweighted events, that helps to keep the storage as small as possible. With the goal of optimising the performance of the unweighting procedure, the so-called partial unweighting is implemented in the Helac-Nlo software, see e.g. [28] for more details. Ntuples allow us to obtain theoretical predictions for different scale choices and PDF sets. Thus, for example the error PDF sets can easily be employed to calculate the internal NLO PDF uncertainties. Furthermore, any infrared-safe (IR-safe) observable can be generated, ranges and bin sizes can be adjusted while no additional time consuming running of the code is required.

LHC setup
We consider the pp → e + ν e µ −ν µ bb bb + X process at NLO in QCD for LHC Run II energy of √ s = 13 TeV. Specifically, α s corrections to the born-level process at O(α 4 s α 4 ) are evaluated. Different lepton generations are used to avoid virtual photon singularities stemming from γ → + − . We do not consider τ leptons as they are usually studied separately at the LHC due to their very rich and complex decay pattern. We keep the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix diagonal. We take LO and NLO NNPDF3.1 PDF sets [55] as the default PDF sets. In both cases they are obtained with α s (m Z ) = 0.118. However, we will show results for other PDF sets, specifically for CT18 [56] and MMHT14 [57]. The difference between various PDF sets originate, among others, from the choice of the data used and the theoretical assumptions made for the global fit. Consequently, it is desirable to see theoretical predictions also for different PDF sets. The running of the strong coupling constant α s with two-loop (one-loop) accuracy at NLO (LO) is provided by the LHAPDF interface [58] with N F = 5. The SM input parameters for our calculations are given by  All other partons, including bottom quarks, and leptons are treated as massless particles. The top quark width is treated as a fixed parameter throughout this work and its value is evaluated for α s (µ R = m t ).
The top quark width at LO has been computed based on the formulas from Ref. [59], whereas the NLO QCD value has been obtained upon applying the relative QCD corrections given in Ref. [60] to the LO width. For our set of input parameters we have Since we treat bottom quarks as massless partons there are no diagrams with Higgs boson exchange at tree level. The Higgs boson contribution appears only at NLO through closed fermion loops involving top quarks. We checked that this contribution amounts to 0.3% of the total NLO cross section for our setup and is therefore negligible. Given its very small impact we decided to neglect the Higgs-boson contribution, which is equivalent to taking the m H → ∞ limit. The electromagnetic coupling α is calculated from the Fermi constant G µ (G µ −scheme) via where G µ is measured in the muon decay. In the G µ −scheme electroweak corrections related to the running of α are taken into account. The evaluation of the residual theoretical uncertainties due to the not yet calculated higher order corrections is based on the exploration of the cross section dependence on the renormalisation scale, µ R , and on the factorisation scale, µ F . For a given definition of the µ 0 scale, judiciously chosen to absorb the large logarithmic corrections that appear at higher orders, we set µ R = µ F = µ 0 but vary the two scales independently in the range with the following additional condition This means that none of the ratios µ F /µ 0 , µ R /µ 0 and µ F /µ R can be larger than two or smaller than one-half. In this way we avoid having in the perturbative expansion logarithms of arguments larger than a chosen amount, regardless of its arbitrariness. In practice, it comes down to considering the following pairs By searching for the minimum and maximum of the resulting cross sections one obtains an uncertainty band. The narrower the band is, the smaller the higher order corrections are expected to be. The variation of the cross section with respect to the scale choice is unphysical. It is just a reflexion of the truncation of the perturbative series. Indeed, if the cross sections are known to all orders, they will not exhibit this dependence. The scale variation is, thus, by no means a rigorous way to estimate the theoretical uncertainty. At best, it might only give an indication of the true uncertainty. For µ 0 we consider two cases. The mass of the heaviest particle appearing in the process is typically considered to be a natural scale choice. Thus, for our fixed scale setting we choose µ 0 = m t . Integrated fiducial cross sections are mostly influenced by final state production relatively close to the tt threshold, which justifies our choice. However, differential cross sections extend up to energy scales that are much larger than the threshold. To this end we define the dynamical scale µ 0 = H T /3. The latter is given by where b i , i = 1, . . . , 4 stands for the four b-jets and p miss T is the missing transverse momentum built out of the two neutrinos (ν e ,ν µ ). This particular choice is especially suitable for the calculations with off-shell top-quark contributions as it carries no information about the underlying top-quark resonant history. Thus, the top-quark reconstruction is not attempted here. The PDF uncertainties together with the scale dependence are the two most important ingredients of the theoretical error on the predictions of the cross sections. All final-state partons with pseudo-rapidity |η| < 5 are recombined into jets with separation R = 0.4 in the rapidity-azimuthal-angle plane via the IR-safe anti−k T jet algorithm [61]. Furthermore, we use the four-momentum recombination scheme. In the first part of the paper contributions induced by the bottom quarks are neglected and we require exactly four b-jets as well as two charged leptons. In the second part of the paper, when the bottom quarks are going to be included in the initial state, we will replace the requirement of having exactly four b-jets with another one, namely, we will ask for at least four b-jets in the final state. We put no restriction on the kinematics of the extra light jet (if resolved) and on the missing transverse momentum. All final state particles and b-jets have to fullfil the following criteria, which we consider to be very inclusive selection cuts

Integrated fiducial cross sections
With the input parameters and cuts specified in Section 3, we arrive at the following predictions for the fixed scale choice µ R = µ F = µ 0 = m t and the default PDF sets At the central scale µ 0 = m t , the gg channel dominates the total LO pp cross section by 94%, followed by the qq+qq channels with 6%. The full pp cross section receives positive and large NLO QCD corrections of 89%. The theoretical uncertainties resulting from scale variation taken in a very conservative way as a maximum of the lower and upper bounds are 65% at LO and 22% at NLO. Therefore, by going from LO to NLO we have reduced the theoretical error by a factor of 3. In the case of truly asymmetric uncertainties sometimes it is more appropriate to symmetrise the errors. After symmetrisation the scale uncertainty at LO is 51% and at NLO does not change substantially, i.e. it is reduced down to 20%. The K-factor that we have obtained K = 1.89, is defined as the ratio of NLO to LO cross sections. In our case both LO and NLO integrated fiducial cross sections are calculated for LO and NLO PDF sets as obtained with α s (m Z ) = 0.118. Had we used the LO NNPDF3.1 PDF set with α s (m Z ) = 0.130 our LO prediction would rather be This would result in the K-factor of K = 1.45. Both findings are correct and reflect the different dependence on the scale of LO and NLO cross sections. Indeed, the LO cross section is much more sensitive to the variation of scales and can change more rapidly than the NLO one. Independently, we can conclude at this point that NLO QCD corrections are large and indispensable to correctly describe the process at hand. Another source of the theoretical error comes from the parameterization of the NNPDF3.1 PDF set. These uncertainties are due to experimental errors in the various data that are used in the fits. They do not, however, take into account additional systematics coming from the underlying assumptions that enter the parametrisation of different PDF sets. The latter cannot simply be quantified within a given scheme. Therefore, we also provide NLO QCD results for two other PDF sets, namely CT18NLO and MMHT2014. We use the corresponding prescription from each group to provide the 68% confidence level (C.L.) PDF uncertainties 2 . Both CT18NLO PDFs and MMHT2014 PDFs include a central set as well as respectively N = 58 and N = 50 error sets in the Hessian representation. The NNPDF3.1 PDF set uses the MC sampling method in conjunction with neural networks. In that case the PDF uncertainties are obtained using the replicas method with a set of N = 100 MC PDF members. The internal NNPDF3.1 PDF uncertainties are very small, at the level of 1% only. Our findings for CT18NLO and MMHT2014 PDF sets are, on the other hand, given by We can see that for CT18NLO and MMHT2014 the internal PDF uncertainties are slightly larger, i.e. of the order of 3%. These uncertainties are nevertheless similar in size to the differences between NLO QCD results obtained with various PDF sets. Indeed, the relative difference between CT18NLO and NNPDF3.1 is 3%, whereas for MMHT2014 and NNPDF3.1 we have instead 1%. We additionally present LO predictions for various LO PDF sets. Due to the lack of LO CT18 PDF sets we went back to the older version and used instead the two LO CT14 PDF sets [62]. More specifically, we employ CT14lo, CT14llo and LO MMHT2014 with α s (m Z ) = 0.118, 0.130 and 0.135 respectively. The larger the value of α s (m Z ) the larger the resulting LO cross section would be. The latter goes to the denominator of the K-factor. Our findings confirm this pattern and can be summarised as follows  In turn, we receive the following spread of the K-factor values 1.81, 1.37, 1.23 respectively. Depending on the LO PDF set employed the NLO QCD corrections to pp → ttbb production in the di-lepton top quark decay channel range from 89% to 23%. On the other hand, the NLO theoretical error is completely dominated by the scale dependence and is consistently at the 22% level (at the 20% level after symmetrisation) independently of the PDF set used. For the dynamical scale setting µ R = µ F = µ 0 = H T /3 our results read (4.5) Our NLO QCD findings for CT18NLO and MMHT2014 PDF sets are given by  As expected, theoretical predictions for µ 0 = H T /3 for the integrated fiducial cross sections are very similar to the results for µ 0 = m t . We show these results for completeness and because we will employ the former scale choice at the differential level, where it plays a crucial role. In Figure 4, the graphical presentation of the behaviour of LO and NLO cross sections upon varying the central value of µ R and µ F by a factor of ξ ∈ {0.125, ..., 8} is shown for the NNPDF3.1 PDF sets. Both cases µ 0 = m t and µ 0 = H T /3 are depicted, that allowed us to compare the two scales. For the sake of completeness, in Figure 4 we present again the scale dependence of the LO and NLO integrated cross sections for each case of µ 0 separately. Also shown is the variation of µ R with the fixed value of µ F and the variation of µ F with fixed µ R . We can note that from the point of view of the integrated fiducial cross sections each scale is a valid choice that might be used in phenomenological applications. We can also add that in the range ξ ∈ {0.5, 1, 2} the scale variation is driven by the changes in µ R . In other words, should we vary µ R and µ F up and down by a factor of 2 around µ 0 simultaneously instead of independently the scale dependence uncertainties would not change significantly. Furthermore, not only the choice of µ 0 but also the value of ξ is important. The latter if not properly selected can introduce large higher order effects and even result in negative (unphysical) NLO cross sections as can be observed in Figure 4 for ξ 0.25. Consequently, such small values should not be selected to set the scale for this process.

Stability test of NLO fiducial cross sections
In Table 3 we show the integrated fiducial cross sections at LO and NLO for different cuts on the transverse momentum of the b-jet. Theoretical uncertainties coming from scale variation are denoted as δ scale , whereas the internal NNPDF3.1 PDF uncertainties are labeled as δ PDF . Also shown is the Kfactor. Results are evaluated using µ R = µ F = µ 0 with µ 0 = m t and for the LO and NLO NNPDF3.1 PDF sets. We observe a very stable behaviour of systematics when varying the p T (b) cut, within the 25 − 40 GeV range. Specifically, δ scale is consistently of the order of 20% and the PDF uncertainties are small for each value of the p T (b) cut. The size of NLO QCD corrections is reduced from 89% to 68%. This reduction of 21% is well within the NLO uncertainties for this process. Results are not changed when they are generated for the dynamical scale choice. Theoretical predictions at LO and NLO for µ 0 = H T /3 are also given in Table 3. For the nominal value of the p T (b) cut the K-factor is slightly larger, i.e. K = 1.94. It is reduced down to K = 1.83 for p T (b) > 40 GeV. The 11% difference is again well within δ scale . We can conclude that the perturbative expansion in α s for the process at hand is not spoiled by the appearance of large logarithms, thus, under excellent theoretical control. 1.83 Table 3. LO and NLO integrated fiducial cross sections for the pp → e + ν e µ −ν µ bb bb + X process at the LHC with √ s = 13 TeV. Results are evaluated using µ R = µ F = µ 0 with µ 0 = m t and µ 0 = H T /3. LO and NLO NNPDF3.1 PDF sets are used. We display results for four different values of the p T (b) cut. Also given are the theoretical uncertainties coming from scale variation (δ scale ) and PDFs (δ PDF ). In the last column the K-factor is shown.

Additional cuts and comparison with ATLAS results
In the following we will examine the behaviour of the integrated cross section upon adding additional cuts, i.e. making the available phase space for the 2 → 8 process more exclusive. We will show results for µ R = µ F = µ 0 = m t only. We have checked, however, that the LO and NLO theoretical predictions for µ 0 = H T /3 are very similar. In the first step we increase p T ( ) from 20 GeV to 25 GeV. Furthermore, we introduce an extra cut, i.e. the separation between the b-jet and the charged lepton in the rapidity-azimuthal-angle plane of ∆R( b) > 0.4. The NLO results with these cuts read The corresponding K-factor is K = 1.88. In the next step another cut has been added, i.e. the separation between charged leptons of ∆R( ) > 0.4. We report the following NLO cross section for this case σ NLO pp→e + νeµ −ν µbbbb (NNPDF3.1, µ 0 = m t ) = 9.52  These two additional conditions for ∆R( ) and p miss T do not affect the K-factor. It remains at the K = 1.88 level. We also notice that the theoretical uncertainties due to scale dependence are remarkably stable for all three cases. A few comments are in order. The first selection that we have applied, namely with the anti − k T jet algorithm and the radius parameter R = 0.4, is very close to the phase-space volume in which the fiducial ttbb cross section has recently been measured by the ATLAS collaboration in the eµ top-quark decay channel [23]. In that paper, among others, σ ttbb was determined by requiring exactly one electron and one muon (with opposite charges) and at least four b-jets. We shall label this final state as eµ + 4b in the following. The measured cross section, after the estimated contributions from ttH, ttW ± and ttZ productions were subtracted, was compared with the theoretical predictions for the ttbb process. The experimental result was found to be higher than the SM prediction but still compatible within the quoted uncertainties. In detail, the measured ATLAS result is given by where 6.5 fb (26%) is the total uncertainty for this measurement, i.e. statistical plus systematical uncertainty. In Table 4 we show various theoretical predictions for ttbb production that have been used in Ref. [23] to compare with the ATLAS data. In particular, results shown there have been generated with the following MC frameworks: Sherpa+OpenLoops [18], PowHel+Pythia 8 [17,19,20] and MadGraph5 − aMCNlo+Powheg-Box+Pythia 8 [43], that are commonly used by the ATLAS and CMS experimental collaborations. In all cases massive b-quarks and four-flavour PDF sets have been employed. The second PowHel result has been obtained for massless b-quarks with the five-flavour PDF set. All theoretical results are based on the NLO matrix element calculations for the on-shell ttbb process. Let us note here, that the measured inclusive fiducial cross section exceed theoretical predictions for the ttbb process. However, they are all compatible within the given uncertainties. For comparison, we also add our result multiplied by a factor of 2 to account for the two decay channels e + µ − and e − µ + . We can observe that the σ result is closer to the σ ATLAS eµ+4b one than other theoretical predictions given in Table 4. Both results agree very well within the quoted uncertainties. Specifically, the agreement of 0.7σ has been obtained 3 . A dedicated comparison between Helac-Nlo and predictions obtained with the help of various NLO matrix element calculations matched to parton shower programs would be in order to understand the source of the spread among these theoretical results. We leave such studies for the near future. We note that in the experimental result as well as in the theoretical predictions used in Ref. [23] also leptonic τ ± decays were included, i.e. e ± and µ ± from τ ± decays were incorporated into the analysis. This is not the case for the Helac-Nlo simulations. Nonetheless, we can roughly estimate the size of the missing contribution by using the NLO fiducial cross section for pp → τ + ν τ τ −ν τ bb bb + X multiplied by the corresponding branching ratio for the leptonic τ + τ − decays. Since our selection cuts are very inclusive this estimate should not be very far away from the true result. Assuming the following branching ratios BR(τ − → µ −ν µ ν τ ) = 17.39% and BR(τ − → e −ν e ν τ ) = 17.82% [63] we can estimate that the Helac-Nlo result should be increased by 3 We note that the scale choice used in Ref. [23] is different than our dynamical scale setting µR = µF = µ0 = HT /3.  which is only 0.6σ away from the ATLAS measured cross section. We conclude this part by saying that the precision of the measurement is still slightly lower than that of the theoretical prediction as obtained with the help of the Helac-Nlo MC program. In the former case the experimental total error is at the 26% level. It comprises statistical and systematical uncertainties. In the latter case the estimated theoretical error is 22%. It combines scale dependence and PDF uncertainties.

Differential fiducial cross sections and PDF uncertainties
In addition to the normalization of the integrated fiducial cross section higher-order QCD corrections can affect the shape of various kinematic distributions. To quantify the size of these effects and to investigate shape distortions that are introduced in this way we shall test a variety of differential distributions for the pp → e + ν e µ −ν µ bb bb + X process. All observables are obtained for the µ R = µ F = µ 0 = H T /3 scale choice, NNPDF3.1 PDF sets and default cuts and parameters. Even though we present results for the kinematic-dependent choice of µ R and µ F only we have also studied differential predictions for the kinematic-independent scale choice, i.e. for µ R = µ F = µ 0 = m t . We will use the latter findings to provide relevant comments. However, in order not to lengthen the manuscript unnecessarily, they will not be shown. This is additionally justified by the fact that at the differential level µ 0 = m t cannot describe efficiently the multi-scale kinematics of the ttbb process. The fixed scale setting leads to perturbative instabilities in the TeV region, where large negative corrections are visible in the tails of several distributions, see e.g. Ref. [28] for a detailed discussion of this behaviour. We show in the upper plots the absolute LO and NLO QCD predictions for the pp → e + ν e µ −ν µ bb bb+X process. The lower panels display the differential K-factors defined as K = dσ NLO (µ 0 )/dσ LO (µ 0 ). Additionally, we provide the uncertainty bands from scale variation defined according to K(µ = ξµ 0 ) = dσ NLO (µ = ξµ 0 )/dσ LO (µ 0 ). The LO blue bands are given to illustrate the relative scale uncertainty of the LO cross section. They are defined according to dσ LO (µ = ξµ 0 )/dσ LO (µ 0 ). The LO and NLO uncertainty bands are obtained by performing a 7-point scale variation around the central value µ 0 .
We note here that our analysis is not trying to identify the origin of b-jets, i.e. it does not distinguish between the extra (prompt) b-jets and b-jets that come from the top-quark decays. This is to avoid the use of the reconstruction techniques for assigning b-jets to a specific production process. We leave such topic for future studies in which modeling uncertainties will be explored in more detail.
We start with the transverse momentum of the four b-jets. They are depicted in Figure 5. The heavy-flavour jets are ordered in p T , thus, p T (b 1 ) corresponds to the leading (hardest) b-jet while p T (b 4 ) is the fourth-leading (softest) b-jet. The differential cross sections are plagued by the same large higher order QCD effects as the integrated fiducial cross sections. Indeed, for example for p T (b 1 ) we observe 90% − 135% NLO QCD corrections, that introduce shape distortions of the order of 45%. The theoretical uncertainties due to scale dependence are, on the other hand, in the range of 20% − 30%. For p T (b 2 ), p T (b 3 ) and p T (b 4 ) higher order QCD corrections are of the order of 70% − 120% while uncertainties are between 10% and 25% depending on the p T of the b-jet. Had we used µ 0 = m t instead we would rather obtain shape distortions up to even 150%. This shows the advantage of the dynamical scale over the fixed one.
In Figure 6 we show the differential cross section distribution as a function of the angular separation between the b-jets, ∆R(bb). Also presented are the transverse momentum and invariant mass of the bb system. They are labelled as p T (bb) and M (bb) respectively. Specifically, we display the two hardest b-jets, denoted as b 1 b 2 , and the two softest b-jets, denoted as b 3 b 4 . Looking at ∆R(b 1 b 2 ) we can notice that the b 1 b 2 system originates predominately from top-quark pair production as b 1 and b 2 are generated mostly in back-to-back configurations. This is additionally confirmed by the p T (b 1 b 2 ) and M (b 1 b 2 ) distributions, that have harder spectra in comparison to the b 3 b 4 system. The latter system is expected to receive large contributions from gluon splittings as manifested by the enhancement at the beginning of the ∆R(b 3 b 4 ) distribution. However, we can also notice rather large contributions in the back-to-back configurations for ∆R(b 3 b 4 ). This suggest that the simple picture that the two high p T b-jets are from top-quark decays while the two low p T b-jets, which are closest in ∆R(bb), are b-jets from the g → bb splitting may not apply. The reconstruction of the production mechanisms for all final states is rather cumbersome when multiple b-jets are present. It requires good reconstruction techniques and excellent understanding of the modelling of top-quark decays. The presence of the additional contributions either from off-shell top quarks or from additional resolved light and/or b-jets makes this picture even more complicated. As we have mentioned earlier we leave such studies on the identification of the origin of the b-jets for the future. Instead, in the following we will focus on the size of NLO QCD corrections to various observables constructed for the b 1 b 2 and b 3 b 4 system. We underline here the fact that such spectra are measured experimentally, see e.g. Ref. [23]. When examining dimensionful observables we notice that the b 1 b 2 system receives larger corrections up to even 150%. For b 3 b 4 , on the other hand, corrections are more constant and between 80% and 110%. For dimensionless angular distributions NLO QCD effects are within the 80%−100% range. All observables, except for p T (b 1 b 2 ), have moderate shape distortions. For the p T (b 1 b 2 ) distribution, however, they are even up to 80%. Had we used the fixed scale choice also in this case we would rather obtain more than 100% changes in the shape of the observables due to QCD higher order corrections. Furthermore, for all observables moderate theoretical uncertainties up to 25% − 30% are observed.
An interesting set of observables is depicted in Figure 7. They are the scalar sums built out of the transverse momenta of the various final states from the pp → e + ν e µ −ν µ bb bb + X process. In particular, we present H T , H vis T , H had T and H lep T . The first one has already been defined as where 1,2 = e + , µ − and p miss T is the missing transverse momentum from the two neutrinos. We also have the visible, hadronic and leptonic versions of H T . The three observables are measured experimentally and defined as follows Although the main emphasis is on the understanding of the b-jet kinematics, we can calculate higher order QCD corrections to any IR-safe observable, which can be constructed from the available final states, in the pp → e + ν e µ −ν µ bb bb + X process. Therefore, we can examine in detail for example the two charged leptons. The advantage of the leptonic observables over hadronic ones lies in the fact that measurements of lepton observables are particularly precise at the LHC due to the excellent lepton energy resolution of the ATLAS and CMS detectors. In Figure 8 we present the following leptonic observables: the transverse momentum of the muon (p T (µ − )), the two leptons' invariant mass (M (e + µ − )), angular difference in the transverse plane (∆φ(e + µ − )), and angular separation (∆R(e + µ − )). NLO QCD corrections to leptonic observables are also substantial in the range of 80% − 150%. Theoretical uncertainties, on the other hand, are again up to 20% − 25% only.
As already pointed out the theoretical uncertainty related to the lack of next-to-next-to-leadingorder QCD corrections is only a fraction of the overall theoretical systematics. To complete our analysis we study differential cross section uncertainties due to the PDF parameterisation including the latest fits by several groups, i.e. NNPDF3.1, MMHT2014 and CT18NLO. We have comparatively examined the impact of PDFs and scale variations on the overall theoretical uncertainty for differential cross section distributions for the pp → e + ν e µ −ν µ bb bb + X process. For all observables that we have examined we could confirm that the PDF uncertainties are small and well below the theoretical uncertainties predicted by scale variation. As an example in Figure 9 we present NLO differential cross section distributions as a function of the transverse momentum of the 1 st , 2 nd , 3 rd and 4 th hardest b-jet. Upper panels show the absolute NLO predictions for three different PDF sets as obtained with the dynamical scale setting. Middle panels display the relative scale uncertainties of NLO predictions for the NNPDF3.1 PDF set. Also presented are the relative predictions for the MMHT2014 and CT18NLO PDF set. They are obtained by normalisation to the central NLO prediction with the NNPDF3.1 PDF set. Finally, lower panels display the relative internal PDF uncertainties of the NLO cross section for each of the three PDF sets. Also here all NLO predictions are normalised to the central NLO prediction with the NNPDF3.1 PDF set. A few comments are in order. To start with, at the differential level the internal PDF uncertainties are up to 11% for CT18NLO and of the order of 8% for MMHT2014. For the default NNPDF3.1 PDF set they are in the 1%−7% range. When comparing the NLO QCD results as obtained with either CT18NLO or MMHT2014 to the findings generated for the default NNPDF3.1 PDF set the relative differences are up to 8% and 3% respectively. Thus, also at the differential level the internal NNPDF3.1 PDF uncertainties are similar in size to differences between the PDF sets that are recommended for applications at the LHC Run II [64].
To summarise this part, for the pp → e + ν e µ −ν µ bb bb + X production process at the LHC for a centre-of-mass-system energy of √ s = 13 TeV and with rather inclusive selection cuts not only big NLO QCD corrections but also significant shape changes are visible when going from LO to NLO. This confirms that NLO QCD effects to this process are extremely important. The theoretical uncertainties due to scale dependence for µ 0 = H T /3 are rather moderate of the order of 20% − 30%. For the fixed scale setting they are much higher. The uncertainties due to the NNPDF3.1 PDF parameterisation are small, i.e. in the 1% − 7% range. Overall, when various PDF sets are examined the PDF uncertainties are maximally up to 11%. Consequently, the dominant component of the final theoretical error for the pp → e + ν e µ −ν µ bb bb + X process is determined by the scale dependence.    Results are shown for the pp → e + ν e µ −ν µ bb bb + X process at the LHC with √ s = 13 TeV. Three different PDF sets are employed. Lower panels display the theoretical uncertainties from scales and PDFs. All predictions are normalised to the central NLO prediction with the NNPDF3.1 PDF set.

Contribution of initial state bottom quarks
In the next step we study the contributions that are induced by the initial state bottom quarks. To this end, additional subprocesses are included in the calculation. For the LO part we add bb → e + ν e µ −ν µ bb bb, bb → e + ν e µ −ν µ bb bb andbb → e + ν e µ −ν µ bbbb. The three reactions are related by crossing symmetry, thus, each comprises 2790 Feynman diagrams. The last two subprocesses must be taken into account as they might not necessarily be suppressed other than by PDFs. Indeed, they are already part of the double-resonant contribution to the tt process with additional two b-jets. For the real emission part, on the other hand, the following subprocesses are included bb → e + ν e µ −ν µ bb bb g, gb → e + ν e µ −ν µ bb bb b, gb → e + ν e µ −ν µ bb bbb, bb → e + ν e µ −ν µ bb bb g andbb → e + ν e µ −ν µ bbbb g. Each subprocess comprises 28728 Feynman diagrams. We continue to use the anti − k T jet algorithm. However a small modification is needed for the jet flavour assignment. As we are dealing with massless bottom quarks from the theoretical point of view the important parton recombination rules, that are required to guarantee the IR-safety of the jet algorithm, are bg → b,bg →b and bb → g. We need, however, additional recombination rules for bb andbb. We employ two variants that are IR-safe at NLO. Beyond NLO IR-safety requires the algorithms of Ref. [65,66].

Charge-blind b-jet tagging
In the first case we use the charge-blind b-jet tagging. From the experimental point of view b-jet tagging algorithms are sensitive mostly to the absolute flavour and they do not additionally tag the charge of the b-jet. In the absence of charge tagging any combination that contains an even number of b and/or b quarks should also be considered to carry zero flavour as from the experimental point of view such signatures will not be distinguishable from bb → g. In this case the complete set of recombination rules for heavy-flavour jets is given by bg → b ,bg →b , bb → g , bb → g ,bb → g . (6.1) We ask for at least four b-jets in the final state and check whether there is at least one combination of (any) four b-jets that fulfils the required cuts. We shall refer to this approach as the charge-blind b-jet tagging or in short as the charge-blind scheme in the following.
From the experimental point of view the disadvantage of this approach might lie in the possibility of the reduction in the b-jet tagging efficiency and in smaller event statistics. In this case the following recombinations rules are employed: We ask for at least four b-jets in the final state, however, this time we check whether any combination with zero total charge, i.e. any bbbb combination, passes the cuts. In this scheme there is no need to consider the bb andbb initiated subprocess. We shall refer to this approach as the charge-aware b-jet tagging or in short as the charge-aware scheme.
In the following the difference between the two approaches is examined. We investigate their impact on the full pp cross sections at LO and NLO in QCD. Specifically, we compare the two schemes to the case when the initial state bottom-quark contribution is omitted. At LO the two approaches are equivalent because we always have exactly the four b-jets that need to pass the cuts. Furthermore, at µ 0 = m t 6.561(2) 0.4367(1) 0.008607(7) 0.006184 (8) µ 0 = H T /3 6.404(3) 0.4092(1) 0.008428(7) 0.006005(7) Table 5. LO integrated fiducial cross section for the pp → e + ν e µ −ν µ bb bb + X production process at the LHC with √ s = 13 TeV for µ 0 = m t and µ 0 = H T /3. Theoretical results with and without the initial state b contributions are provided for the LO NNPDF3.1 PDF set. Also given are Monte Carlo integration errors (in parenthesis).  Table 6. NLO integrated fiducial cross section for the pp → e + ν e µ −ν µ bb bb + X production process at the LHC with √ s = 13 TeV for µ 0 = m t and µ 0 = H T /3. Theoretical results with and without the initial state b contributions are provided for the NLO NNPDF3.1 PDF set. Results for four different values of the p T (b) cut and for a jet veto with p veto T (j) = 50, 100 GeV are additionally provided. Also given are the theoretical uncertainties due to scale variation (δ scale ) and PDFs (δ PDF ). Finally, Monte Carlo integration errors are shown (in parenthesis).
this order both schemes differ only by the initial state subprocesses that must be taken into account. Simply, in the charge-aware scheme the bb andbb subprocesses are not considered.
In Table 5  In Table 6  A few comments are in order. If we look at the bottom-quark initiated subprocesses only, their combined contribution at NLO is dominated by the bg/bg initial states. Furthermore, at this level there are substantial differences between the two schemes. Specifically, for µ 0 = m t we have for the charge-aware and charge-blind scheme the following results where in bracket the MC error is additionally provided. It is only when the dominant gg partonic subprocess is taken into account that the differences between the two schemes become insignificant. Furthermore, having at hand the results with an additional jet veto on the fifth jet we can make the following remarks on the size of the higher order QCD effects. We observe that the jet veto can reduce the K-factor for the default set of cuts. Specifically, for µ 0 = m t (µ 0 = H T /3) when p veto T (j) = 100 GeV is used we obtain K = 1.48 (K = 1.58), whereas for p veto T (j) = 50 GeV we have K = 1.11 (K = 1.23). As a bonus of our study in Figure 10 we show graphically NLO integrated fiducial cross section for the pp → ttbb process with and without the initial state b contributions. Results are given for the NLO NNPDF3.1 PDF set. However, we also provide results with other PDF sets for the case when the initial state b-quark contributions are not included. We display the result for CT14, CT18, NNPDF3.0, MMHT2014 and for the NLO ABMP16 PDF set [71]. In the latter case the internal PDF uncertainties are evaluated using the symmhessian method, see e.g. Ref. [58]. To put all the results in perspective the theoretical uncertainty resulting from scale dependence are also shown. All predictions are calculated for µ 0 = H T /3. We notice that, all theoretical predictions agree within their internal  Finally, we perform a similar comparative analysis at the differential level. In Figure 11 we display NLO differential cross section distributions as a function of p T (b 1 ), ∆R(b 1 b 2 ), M (b 1 b 2 ) and p T (b 1 b 2 ). The upper plots show absolute NLO QCD predictions with and without the initial state b-quark contributions. Also shown is the NLO scale dependence for the default case. The lower panels display the differential ratios of the results with the initial state bottom-quark contributions to the one without such contributions. Overall, for the process at hand when comparing to the differential theoretical errors due to scale dependence the inclusion of the bottom-quark induced subprocesses is not important also at the differential level. Nonetheless, a consistent treatment of heavy-flavour jets is necessary to obtain IR-finite results.

Comparison with previous results
In what follows we compare our predictions for the pp → e + ν e µ −ν µ bb bb + X process at the LHC running at √ s = 13 TeV to the theoretical predictions from Ref. [25], which we dub as DLP computation. The comparison is carried out at the integrated and differential level at LO and NLO in QCD. We adopt the dynamical scale setting used in Ref. [25]. Specifically, we employ µ R = µ F = µ 0 where µ 0 = µ DLP is given by Furthermore, E T (i) = p 2 T (i) + M 2 (i) and M 2 (i) is the invariant mass squared of the object considered. Similarly to our dynamical scale setting also in this case the top-quark reconstruction is not for the pp → e + ν e µ −ν µ bb bb + X production process at the LHC with √ s = 13 TeV. The upper plots show absolute NLO QCD predictions without and with the initial state bottom quark contributions. In the latter case results are shown for the charge blind and charge aware b-jet tagging. The lower panels display the differential ratios of these contributions. Also shown is the NLO scale dependence for the case without the initial state bottom quark contributions. The NLO NNPDF3.1 PDF set is employed and µ 0 = H T /3 is used.
Scale dependence of the LO and NLO integrated fiducial cross sections for the pp → e + ν e µ −ν µ bb bb + X production process at the LHC with √ s = 13 TeV for µ 0 = µ DLP . The LO and NLO NNPDF3.1 PDF sets are employed. Also shown is the variation of µ R with fixed µ F and the variation of µ F with fixed µ R as well as the comparison to our default scale choice µ 0 = H T /3.
It can be directly compared to the result form Ref. [25] σ LO DLP (NNPDF3.1, µ 0 = µ DLP ) = 5.198(4) +60% −35% fb . 3) The two results agree perfectly within the given MC errors. Had we additionally included the two missing subprocesses bb andbb the result would rather be σ LO HELAC−NLO = 5.206(2) fb. The latter is 1.8σ away from σ LO DLP . In order to reproduce the NLO QCD result from Ref. [25] the additional fifth jet (if resolved) must be added to the definition of µ 0 = µ DLP . The µ 0 = µ DLP scale is, thus, given by The Helac-Nlo NLO result and the σ NLO DLP one, both without the bb andbb contributions, read  Theoretical predictions are in perfect agreement. Had we also included the bb andbb contributions at NLO our result would rather be σ NLO HELAC−NLO = 10.29(1) fb, so still in agreement with σ NLO DLP . A few comments are in order. Both LO and NLO QCD results for µ 0 = µ DLP are much smaller than these obtained either for µ 0 = m t or µ 0 = H T /3. In the case of the NLO predictions the differences exceed the theoretical uncertainties due to scale dependence. The dynamical scale µ 0 = µ DLP includes the dependence on the p T of the additional jet, thus, on the average its value is larger. On the other hand, the asymptotic freedom guarantees that the value of α s becomes smaller for larger µ 0 , resulting in lower NLO and LO cross sections. Had we removed the additional jet from the definition of µ 0 = µ DLP the resulting NLO cross section would increase by about 13%. Consequently, it would be in agreement with the NLO findings both for µ 0 = m t and µ 0 = H T /3 within the quoted theoretical uncertainties.  Figure 13. Differential cross section distributions as a function of p T (b 2 ), p T (µ − ), ∆φ(e + µ − ) and cos θ(e + µ − ) for the pp → e + ν e µ −ν µ bb bb + X production process at the LHC with √ s = 13 TeV. Results are presented at LO and NLO for µ 0 = µ DLP . Comparison between Helac-Nlo predictions and results from Ref. [25] (given with the subscript DLP) is shown. The NNPDF3.1 PDF sets are employed. Also displayed are Monte Carlo integration errors.
In Figure 12 we provide the graphical representation of the scale dependence of the LO and NLO integrated fiducial cross sections for the pp → e + ν e µ −ν µ bb bb + X production process as obtained with µ 0 = µ DLP . Also shown are the variation of µ R with fixed µ F and the variation of µ F with fixed µ R . Finally, the comparison to our dynamical scale setting of µ R = µ F = µ 0 = H T /3 is depicted.
In Figure 13 we depict the comparison between the Helac-Nlo predictions and the theoretical results from Ref. [25] at the differential level. Specifically, we show LO and NLO differential cross section distributions as a function of p T (b 2 ), p T (µ − ), ∆φ(e + µ − ) and cos θ(e + µ − ). Also here heavyflavour induced subprocesses are included except from bb andbb. For each observable good agreement has been found.

Summary
In this paper we provided state-of-the-art theoretical predictions for the pp → e + ν e µ −ν µ bb bb + X process at the LHC with √ s = 13 TeV. In the computation off-shell top quarks were described by Breit-Wigner distributions. Furthermore double-, single-as well as non-resonant top-quark contributions and interference effects were consistently incorporated at the matrix element level. Non-resonant and off-shell effects due to the finite W -boson width were also consistently taken into account. All results were obtained with the help of the Helac-Nlo MC package. We have shown LO and NLO predictions for the integrated and differential cross sections, that are phenomenologically relevant for LHC physics. We assessed the theoretical uncertainties of our high-precision theoretical predictions stemming from scale dependence and PDF parameterisation while using fixed and dynamical scale settings, i.e. µ R = µ F = µ 0 , where µ 0 = m t or µ 0 = H T /3. Furthermore, the following PDF sets were examined NNPDF3.1, CT18 and MMHT14. The full pp cross section receives positive and large NLO QCD corrections of 89%. The theoretical uncertainties resulting from scale variation are 65% at LO and 22% at NLO. The internal PDF uncertainties are very small, at the level of 1% − 3% only. Thus, the NLO theoretical error is completely dominated by the scale dependence. We provide LO and NLO results for other PDF sets and observe that the K-factor has a very large spread from 1.81 down to 1.23 depending on the LO PDF set employed and specifically on the value of the strong coupling constant α s (m Z ) used. On the other hand, we observed a very stable behaviour of the systematics when varying the p T (b) cut or adding additional cuts.
The differential cross sections have been plagued by the same large higher order QCD effects as the integrated fiducial cross sections. Not only big NLO QCD corrections but also significant shape changes were visible when going from LO to NLO. This confirms that NLO QCD effects to the pp → e + ν e µ −ν µ bb bb + X process are extremely important. The theoretical uncertainties due to scale dependence for µ 0 = H T /3 are rather moderate of the order of 20% − 30%. For the fixed scale setting they are much higher. The uncertainties due to the NNPDF3.1 PDF parameterisation are small, i.e. in the 1% − 7% range. When other PDF sets have been examined the PDF uncertainties increased maximally up to 11%. Consequently, the final theoretical error for the process at hand remains dominated by the scale dependence.
In the next step we have studied the contributions that are induced by the initial state bottom quarks. To this end, additional subprocesses were included in the calculation. Additionally, we needed new recombination rules for partons to construct light-and heavy-flavour jets. We employed two variants that are IR-safe at NLO: charge-blind and charge-aware b-jet tagging. The differences between the two approaches were examined in detail. Overall, for the process at hand when comparing to the theoretical errors due to scale dependence the size of the contributions generated by the bottom-quark induced subprocesses is negligible both at the integrated and differential level. Nonetheless, a consistent treatment of heavy-flavour jets is necessary to obtain IR-finite results.
Finally, we compared our predictions for the pp → e + ν e µ −ν µ bb bb + X process to the previous results obtained in Ref. [25]. After clarifying the scale choice used in Ref. [25] with the authors and the status for the two bb andbb subprocesses perfect agreement has been found both at the integrated and differential level.