New NLOPS predictions for $\boldsymbol{t\bar{t}+b}$-jet production at the LHC

Measurements of $t\bar{t} H$ production in the $H\to b\bar{b}$ channel depend in a critical way on the theoretical uncertainty associated with the irreducible $t\bar{t}+b$-jet background. In this paper, analysing the various topologies that account for $b$-jet production in association with a $t\bar{t}$ pair, we demonstrate that the process at hand is largely driven by final-state $g\to b\bar{b}$ splittings. We also show that in five-flavour simulations based on $t\bar{t}+$multi-jet merging $b$-jet production is mostly driven by the parton shower, while matrix elements play only a marginal role in the description of $g\to b\bar{b}$ splittings. Based on these observations we advocate the use of NLOPS simulations of $pp\to t\bar{t}b\bar{b}$ in the four-flavour scheme, and we present a new POWHEG generator of this kind. Predictions and uncertainties for $t\bar{t}+b$-jet observables at the 13 TeV LHC are presented both for the case of stable top quarks and with spin-correlated top decays. Besides QCD scale variations we consider also theoretical uncertainties related to the POWHEG matching method and to the parton shower modelling, with emphasis on $g\to b\bar{b}$ splittings. In general, matching and shower uncertainties turn out to be remarkably small. This is confirmed also by a tuned comparison against Sherpa+OpenLoops.


Introduction
At the Large Hadron Collider, searches for ttH production in the H → bb channel are plagued by a large QCD background, which is dominated by ttbb production, and the availability of precise theoretical predictions for this multi-particle background process is of crucial importance for the sensitivity of ttH(bb) analyses. The process pp → ttbb is also very interesting on its own, as it provides a unique laboratory to explore the QCD dynamics of heavy-quark production and to test state-of-the-art Monte Carlo predictions in a nontrivial multi-scale environment.
As a result of its α 4 S dependence, the leading-order (LO) ttbb cross section is highly sensitive to variations of the renormalisation scale. The uncertainty corresponding to standard factor-two scale variations amounts to 70-80% at LO, and the inclusion of next-to-leading order (NLO) QCD corrections [1][2][3] is mandatory. At NLO, the scale dependence goes down to 20-30%, and in order to avoid excessively large K-factors and potentially large corrections beyond NLO, the renormalisation scale should be chosen in a way that accounts for the fact that the typical energies of the b-jet system are far below the hardness of the underlying pp → tt process [3].
The first NLOPS simulation of pp → ttbb was carried out in Powhel [4,5] by combining NLO matrix elements in the five-flavour (5F) scheme with parton showers by means of the Powheg method [6,7]. Shortly after, an NLOPS generator based on four-flavour (4F) pp → ttbb matrix elements became available in the Sherpa+OpenLoops framework [8], which implements an improved version [9] of the MC@NLO matching method [10]. Thanks to the inclusion of b-mass effects, ttbb matrix elements in the 4F scheme are applicable to the full b-quark phase space, including regions where one b-quark remains unresolved. Thus the 4F scheme guarantees a consistent NLOPS description of inclusive tt + b-jet production with one or more b-jets. On the contrary, NLOPS ttbb generators based on 5F matrix elements with massless b-quarks suffer from collinear g → bb singularities that require ad-hoc restrictions of the physical phase space through generation cuts.
In Ref. [8] it was pointed out that matching and shower effects play an unexpectedly important role in tt + b-jet production. This is due to the fact that two hard b-jets can arise from two hard jets involving each a collinear g → bb splitting. In NLOPS simulations of pp → ttbb, such configurations result from the combination of a g → bb splitting that is described at NLO accuracy through ttbb matrix elements together with a second g → bb splitting generated by the parton shower. The impact of this so-called doublesplitting mechanism can have similar magnitude to the ttH(bb) signal, and the thorough understanding of the related matching and shower uncertainties is very important for ttH analyses.
A first assessment of NLOPS uncertainties was presented in Ref. [11] through a tuned comparison of NLOPS ttbb simulations in Powhel [4,5], Sherpa+OpenLoops [8] and Madgraph5aMC@NLO [12]. On the one hand, this study has revealed significant differences between the two generators based on the MC@NLO matching method 1 , i.e. Sherpa and Madgraph5aMC@NLO. Such differences were found to be related to a pronounced dependence on the shower starting scale in Madgraph5aMC@NLO. On the other hand, in spite of the fact that Sherpa+OpenLoops and Powhel implement different matching methods and different parton showers, the predictions of these two generators turned out to be quite consistent. However, due to the limitations related to the use of the 5F scheme in Powhel-which have been overcome only very recently with the 4F upgrade of Powhel [13]-the agreement between Powhel and Sherpa+OpenLoops did not allow to draw any firm conclusion in the study of Ref. [11].
To date the assessment of theoretical uncertainties in ttbb production remains an important open problem. In this context one should address the question of which of the various NLOPS methods and tools on the market are more or less appropriate to describe the process at hand. Moreover, in order to address such issues in a systematic way, it is 1 Here one should keep in mind that the Madgraph5aMC@NLO and Sherpa implementations are not identical. For instance, only the latter guarantees an exact O(αS) description of soft radiation, which is achieved by upgrading the first shower emission and the related Catani-Seymour counterterm to full-colour accuracy.
desirable to develop a better picture of the QCD dynamics that drive tt + b-jet production. In this spirit, this paper starts with a discussion of the various possible frameworks for theoretical simulations of tt + b-jet production at NLOPS accuracy. In particular, we present detailed studies on the role of g → bb splittings and discuss the advantages and disadvantages of the 4F and 5F schemes. To this end we quantify the relative importance of g → bb splittings of initial-state (IS) and final-state (FS) type by using approximations based on collinear QCD factorisation, as well as by decomposing pp → ttbb matrix elements into diagrams involving IS and FS g → bb splittings. These studies demonstrate that tt + b-jet production is widely dominated by pp → ttg followed by FS g → bb splittings. This holds also for observables where initial-state splittings are expected to be enhanced, such as in regions with a single resolved b-jet. These findings support the use of NLOPS generators based on pp → ttbb matrix elements in the 4F scheme, where b-mass effects guarantee a consistent treatment of FS g → bb splittings. We also consider more inclusive simulations of tt+jets production based on multi-jet merging [14][15][16][17][18][19]. In this case we find that tt + b-jet observables suffer from an unexpectedly strong dependence on the parton-shower modeling of g → bb splittings.
Motivated by the above findings we present a new Powheg generator for pp → ttbb in the 4F scheme. 2 At variance with the Powhel generator of Ref. [13], this new Powheg generator is implemented in the Powheg-Box-Res framework [22] using OpenLoops, which guarantees a very fast evaluation of the required 2 → 4 and 2 → 5 matrix elements. The new generator supports also top-quark decays including spin-correlation effects. Moreover, in order to guarantee a more consistent resummation of QCD radiation, the separation of the so-called singular and finite parts in the Powheg-Box is not restricted to initialstate radiation as in Ref. [13] but is applied also to final-state radiation (see Section 3).
For what concerns the Powheg methodology we pay particular attention to issues related to the multi-scale nature of the process at hand. In particular we point out that the treatment of the recoil associated with NLO radiation can induce sizeable distortions of the underlying ttbb cross section. This technical inconvenience restricts the domain of applicability of QCD factorisation in a way that can jeopardise the efficiency of event generation and can also lead to unphysical resummation effects. Fortunately, such issues can be avoided by means of a Powheg-Box mechanism that restricts the resummation of real radiation to kinematic regions where QCD factorisation is fulfilled within reasonably good accuracy.
Predictions for pp → tt + b-jets at the 13 TeV LHC are presented for various cross sections and distributions with emphasis on the discussion of theoretical uncertainties. Besides QCD scale variations, also uncertainties related to the matching method and intrinsic shower uncertainties are analysed in detail. In particular, we consider different approximations for the modelling of g → bb splitting as well as α S scale uncertainties in Pythia. Moreover we compare Pythia to Herwig. Finally, to gain further insights into the size and the nature of matching and shower uncertainties we present a consistent comparison   Figure 1. Sample ttbb diagrams with IS (left) and FS (right) g → bb splittings. In NLOPS simulations of inclusive tt production in the 5F scheme the black subtopologies are described in terms of tree matrix elements, while the orange lines correspond to parton shower emissions.
of Powheg+Pythia generators of ttbb and inclusive tt production, against corresponding generators based on Sherpa+OpenLoops.
The new ttbb generator will be soon publicly available on the Powheg-Box web site [23].
The paper is organised as follows. In Section 2 we study the role of g → bb splittings in tt + b-jet production, and we point out the advantages of Monte Carlo generators based on pp → ttbb matrix elements in the 4F scheme as compared to more inclusive generators of tt production in the 5F scheme. Technical aspects of the new Powheg generator and the setup for numerical simulations are discussed in Section 3. In particular, in Section 3.1 we review aspects of the Powheg method that can play a critical role for multi-scale process like pp → ttbb. Detailed predictions and uncertainty estimates for cross sections and distributions for pp → ttbb with stable and unstable top quarks can be found in Sections 4 and 5, respectively. Our main findings are summarised in Section 6.
2 Anatomy of tt + b-jet production and g → bb splittings Events with tt + b-jets final states arise from an underlying pp → tt process that takes place at scales of the order of 500 GeV and is accompanied by the production of b-jets at typical transverse momenta of a few tens of GeV. The production of b-jets is governed by IS or FS g → bb splittings and is enhanced in kinematic regions where the p T of individual b-quarks becomes small or the bb pair becomes collinear and, possibly, also soft.
The understanding of the QCD dynamics that governs b-jet production is a crucial prerequisite for a reliable theoretical description of tt + b-jet production and related uncertainties. In this spirit, this section compares various theoretical frameworks for the description of tt + b-jet production at NLO QCD accuracy with a special focus on the role of g → bb splittings. Specifically, we compare inclusive or merged simulations of tt+multijet production in the 5F scheme against a description based on ttbb matrix elements in the 4F scheme, pointing out the advantages of the latter.
Numerical studies presented in this section are based on the setup specified in Sections 3.2 and 3.6 and have been performed with Sherpa+OpenLoops.

NLOPS tt simulations in the five-flavour scheme
Inclusive NLOPS generators of tt production [24,25] are based on pp → tt NLO matrix elements matched to partons showers in the 5F scheme. In this framework, as illustrated in Fig. 1, tt + b-jet events are generated starting from 2 → 3 tree matrix elements of type gb → ttb or gg/qq → ttg. In the latter case, ttbb events arise via FS g → bb shower splittings. Instead, in the case of gb → ttb the final-state b-quark emerges from the matrix element, while g → bb splittings generate the initial-state b-quark through the evolution of the 5F PDFs. The unresolved spectator b-quark associated with such IS g → bb splittings is emitted by the parton shower via backward evolution. The main advantage of the 5F scheme lies in the resummation of potentially large α S ln(m t /m b ) terms associated with the evolution of the b-quark density. However, such logarithmic effects are typically rather mild at the LHC [26]. Moreover, as we will show in Section 2.3, tt + b-jet production is largely dominated by topologies with FS g → bb splittings. For this reason, tt + b-jet predictions based on NLOPS tt generators suffer from the twofold disadvantage given by the direct dependence on the parton-shower modelling of FS g → bb splittings plus the LO nature of the underlying ttg matrix element.

tt+multi-jet merging in the five-flavour scheme
As a possible strategy to reduce the sensitivity to the parton shower and increase the accuracy of theoretical predictions we consider tt+multi-jet merging in the 5F scheme. In this approach, a tower of NLOPS simulations for tt + 0, 1, . . . , N jet production is merged into a single inclusive sample [17][18][19]. This is achieved by clustering QCD partons into jets with a certain k T -resolution, Q cut , which is known as merging scale. At LO, the phase-space regions with N = 0, 1, . . . , N max resolved jets (k T > Q cut ) are described in terms of tt + N -jet LOPS simulations. The LOPS simulation with N max jets fills also the phase space with N > N max resolved jets by means of the parton shower. At NLO, the resolution criterion used to separate regions of different jet multiplicity exactly the same as for LO, while the basic difference with respect to LO merging lies in the fact that tt + N -jet LOPS simulations are replaced by corresponding NLOPS simulations. Thus in NLO (LO) merging the effective number of resolved jets that is described at NLOPS (LOPS) accuracy is N eff = min{N, N max }, while the (N eff + 1) th resolved or unresolved jet is described at LOPS (pure PS) accuracy, and all remaining resolved or unresolved jets are described at pure PS accuracy.
Multi-jet merging for tt+jet at NLO can be performed in a fully automated way within the Sherpa [27] and Madgraph5aMC@NLO [12] frameworks. However, the fact that tt+ b-jet events constitute only a small fraction of a tt+jets sample poses very high requirements in terms of Monte Carlo statistics. Moreover, in order to minimise the dependence on parton-shower modelling, such simulations should be performed using a small merging scale and including a sufficiently high number of NLO jets, N max . The required CPU resources grow very fast at large N max and small Q cut , and state-of-the-art merged simulations can handle up to two jets at NLO [28] Figure 2. Sample diagrams describing the interplay between matrix elements (black) and parton shower (orange) for tt + b-jet events in merged LOPS simulations of tt + 0, 1, 2 jet production in the 5F scheme. In regions where the bb system and its parent gluon are produced below the merging scale, the event is described through hard tt matrix elements plus parton-shower branchings (a). When the parent gluon becomes harder, the parton shower is used only to model collinear g → bb splittings (b). Above the merging scale, if g → bb splittings do not belong to the two hardest branchings they are still left to the parton shower, while 2 → 4 matrix elements are used to account for other light jets (c). Otherwise, the event is described through ttbb matrix elements (d).
The multi-jet merging description of ttbb events with FS g → bb splittings is sketched in Fig. 2 for the case of tt+0, 1, 2-jet merging at LO. In regions where the bb pair and/or the parent gluon are emitted at small scales, g → bb splittings are expected to be generated by the parton shower, while hard b-jet pairs are expected to arise from ttbb matrix elements. However, as we will see, typical tt+b-jet events involve additional light jets that are emitted at harder scales with respect to the g → bb branching. In that case, 2 → 4 matrix elements account only for light jets, and g → bb splittings are left to the parton shower. In general, the relative importance of matrix elements and parton shower depends on the resolution scale Q cut , and using a finite resolution is mandatory in the 5F scheme, since collinear g → bb splittings are divergent, i.e. ttbb matrix elements cannot be used in the full phase space.
In Figures 3-4 we analyse the matrix-element content of a tt + 0, 1, 2-jet merged simulation at LO. These studies are based on the MEPS@LO method [16] implemented in Sherpa, but the main findings are expected to hold also for other merging methods. Besides the b-jet multiplicity distribution, in Figures 3-4 we plot differential observables in the presence of ttbb cuts, i.e. requiring N b ≥ 2 b-jets as defined in Section 3.6. For jets we apply the acceptance cuts (3.29) and, in order to maximise the possibility to resolve jets at matrix-element level, we choose a merging scale lower than the jet-p T threshold, Q cut = 20 GeV.
As expected, in Fig. 3 we find that tt + b-jet observables with N b ≥ 2 resolved bjets are largely dominated by tt+2-parton matrix elements. This holds also for N b ≥ 1. However, the breakdown of the merged sample into contributions from matrix elements with different b-quark multiplicity in Fig. 4 reveals that, in spite of the low merging scale, the cross section for producing one or more b-jets is dominated by matrix elements with zero b-quarks. The contribution of ttbb matrix elements hardly exceeds 50% even in the region of large b-jet p T or large invariant mass of the b-jet pair. This counterintuitive  feature can be attributed to the fact that, in tt+jet events that involve g → bb splittings, the two hardest QCD branchings are typically associated with the emission of the parent gluon of the bb pairs and/or with the production of other light jets. As a consequence, in LOPS merged samples with N max = 2, g → bb splittings are left to the parton shower. We have verified that the contribution of ttbb matrix elements remains relatively low even when the merging procedure is extended up to 3 or 4 jets with Q cut = 20 GeV. Moreover, we have checked that increasing the merging scale leads to a further suppression of the contribution of ttbb matrix elements. Fig. 4d demonstrates that tt + 2 b-jet events are indeed accompanied by abundant emission of extra light jets, and the importance of ttbb matrix elements decreases with increasing light-jets multiplicity. Moreover, even in the bin with zero additional light jets it turns out that the contribution of ttbb matrix elements remains below 50%. This is probably due to the fact that g → bb splittings tend to take place at branching scales  below the jet-p T -threshold of 25 GeV.
In summary, contrary to naive expectations, tt+jets samples based on LOPS merging do not guarantee a matrix-element description of b-jet production but largely rely on the parton shower modelling of g → bb splittings. In the case of multi-jet merging at NLO, to a certain extent g → bb shower splittings should be matched to ttbb and ttbbg tree matrix elements. Nevertheless, based on the above observations, the theoretical accuracy in the description of b-jet production is expected to remain between the LOPS and the pure PS level. 3 Finally, in the light of the presence of abundant light-jet radiation with a typical hardness beyond the one of b-jets, the role of hard radiation on top of the ttbb system should Figure 5. Generic leading-order gg → ttbb topologies. The first line shows the most general form of topologies with IS (a) and FS (b) g → bb splittings. The second line shows the generic form of those topologies with IS (c) and FS (d) splittings that turn out to dominate gg → ttbb. The labels ij = 56, 65 stand for the bb system, while α = 1, 2 indicates the initial-state gluon that generates the radiation.
be studied with great care also in the context of NLOPS simulations of ttbb production (see Section 3.1).

ttbb production in the four-flavour scheme
In order to minimise the dependence on parton-shower modelling and to maximise the use of higher-order matrix elements, in the following we will adopt a description of tt + b-jet production based on ttbb matrix elements in the 4F scheme. In this scheme, b-quarks are treated as massive partons, and g → bb splittings are free from collinear singularities. Thus ttbb matrix elements can be used in the entire phase space. Generic ttbb topologies where b-quarks emerge from IS and FS splitting processes are illustrated in Fig. 5. In the case of FS g → bb splittings, ttbb matrix elements with m b > 0 can be extended to the collinear regime, where the bb pair becomes unresolved within a single b-jet. Similarly, 4F ttbb matrix elements describe also collinear IS g → bb splittings, where the spectator b-quark is emitted in the beam direction and remains unresolved, while the bg → ttb sub-process with a single b-jet corresponds to the description of tt + b-jet production at LO in the 5F scheme. Thus, ttbb matrix elements provide a fully inclusive description of tt + b-jet production, and NLO predictions in the 4F scheme yield NLO accuracy both for observables with two b-jets and for more inclusive observables with a single resolved b-jet. The inclusion of m b effects in g → bb splittings represents a clear advantage of the 4F scheme with respect to the 5F scheme. However, the 4F scheme has the disadvantage that potentially large α S ln(m b /Q) terms that arise from IS g → bb splittings are not resummed through the PDF evolution. In the following, in order to assess the relevance of this limitation, we decompose the pp → ttbb LO cross section into contributions from IS and FS g → bb splittings. Since the qq channel involves only FS g → bb splitting, we focus on the gg channel and we first consider a naive diagrammatic splitting of the gg → ttbb matrix element, The terms M IS,ttbb and M FS,ttbb correspond, respectively, to 18 diagrams with IS g → bb splittings and 16 diagrams with FS g → bb splittings. Generic diagrams with IS and FS splittings are depicted in Fig. 5a-b. The term M rem,ttbb corresponds to two remaining diagrams where an s-channel gluon splits into an on-shell b-quark and an off-shell b-line coupled to the tt system. Its numerical impact turns out to be negligible. Based on (2.1) we split the ttbb cross section into three terms, where the IS and FS parts are defined as while dσ int,ttbb consists of the interference between M IS,ttbb and M FS,ttbb plus a minor contribution from M rem,ttbb . In order to check the soundness of the above gauge-dependent separation we compare it to an alternative definition of IS and FS g → bb contributions based on the collinear limits of the gg → ttbb matrix element. In this case we define Here |M IS⊗tt | 2 and |M FS⊗tt | 2 describe the collinear limits of the topologies depicted in Fig. 5c and 5d, respectively. Note that, for simplicity, we consider only the leading collinear enhancements where the bb system originates either through the combination of g → bb and b → gb IS splittings (M IS⊗tt ) or via IS g → gg plus FS g → bb splittings (M FS⊗tt ). For events with external momenta where α ∈ {1, 2} and ij ∈ {56, 65} specify, respectively, the IS gluon emitter and the ordering of the bb pair as depicted in Fig. 5c and 5d. K IS/FS (p α , p i , p j ) are the corresponding splitting kernels. The choice of α and ij specifies a particular topology, and the maximum in (2.6) defines M IS⊗tt 2 and M FS⊗tt 2 as the collinear limit of the most likely topology of IS and FS type. The splitting kernels read with T R = 1/2, C F = 4/3 and C A = 3. The various momentum fractions are set to Finally, the underlying gg → tt squared matrix element in (2.6) reads where helicity/colour sums and average factors are included, and the momentum of the IS emitter has to be rescaled by z.
Numerical results for the diagrammatic decomposition (2.2) and the collinear decomposition (2.4) of pp → ttbb at √ s = 13 TeV are shown in Fig. 6. The first two plots display the leading b-jet p T distribution in the presence of ttb and ttbb cuts as defined in Section 3.6, i.e. requiring N b ≥ 1 and N b ≥ 2 b-jets, respectively. In the ttbb phase space, topologies with g → bb FS splittings turn out to be surprisingly close to the full matrix element, with deviations that do not exceed the 10% level in the entire spectrum. This agreement remains remarkably good also in the inclusive ttb phase space, where IS splitting processes with one unresolved b-quark are expected to be more pronounced. Actually, with ttb and ttbb cuts the pure IS contribution ranges between 20-40% and 15-75%, respectively, but is almost entirely cancelled by the IS-FS interference.
The fact that the collinear approximations (2.4) agree rather well with the corresponding squared Feynman diagrams up to relatively high p T confirms that ttbb production is dominated by the topologies in Fig. 5c and 5d. On the other hand, the importance of interference effects provides strong motivation for using exact ttbb matrix elements, while collinear approximations such as those in (2.6) or in the parton-shower modelling of g → bb splittings should be used with due caution.  Figure 6. Breakdown of the pp → ttbb cross section into contributions from topologies with IS and FS splittings at fixed-order LO in the 4F scheme: distributions in the p T of the leading b-jet with ttb (a) and ttbb (b) cuts, and distributions in the invariant mass (c) and ∆R separation (d) of the two leading b-jets with ttbb cuts. The complete gg/qq → ttbb matrix-element prediction (solid red) is split according to (2.2) into contributions from topologies of IS (solid blue) and FS (solid green) type and their interference (solid purple). This is compared to the gauge-invariant breakdown (2.4) into IS (dashed blue) and FS (dashed green) parts based on the collinear limits of the ttbb matrix element. Note that the qq channel consists solely of FS g → bb contributions.
The above considerations apply also for the m bb and ∆R bb distribution in Fig. 6c and d. In particular, we observe that topologies with FS g → bb splittings are very close to the full matrix element in the whole m bb spectrum as well as for ∆R bb < 2 . At the same time, for 50 GeV < m bb < 200 GeV and 1 < ∆R bb < 2.5, i.e. in the range of interest for ttH(bb) analyses, we observe that IS splitting contributions and negative interference effects grow fast and tend to become very sizable. Thus a naive separation into contributions from IS and FS splittings is not applicable at large m bb and ∆R bb . On the other hand, in the region of moderate invariant mass and ∆R separation, which contains the bulk of the ttbb cross section, interference effects are rather small, and bb pairs turn out to originate almost entirely from FS g → bb splittings.
In summary, given that the 5F scheme is based on the LO process gb → ttb, where FS g → bb splittings and interference effects are entirely neglected, the above observations provide strong motivation for a description of tt + b-jet production based on ttbb matrix elements in the 4F scheme.

Technical aspects and setup of NLOPS ttbb simulations
In this section we introduce a new Powheg generator based on ttbb matrix elements in the 4F scheme. Special emphasis is devoted to some technical aspects of the Powheg method that turn out to play an important role for multi-scale processes like pp → ttbb. In addition we describe the setup used for the tt + b-jet simulations presented in Sections 4-5, i.e. all relevant input parameters, scale choices and parton shower settings, as well as the treatment of theoretical uncertainties and the definitions of physics objects and selection cuts. Finally we provide details on the treatment of top-quark decays.
The new ttbb generator is implemented in the Powheg-Box-Res framework [22], and the relevant LO and NLO matrix elements are computed by OpenLoops [29][30][31] through its Powheg-Box-Res interface [32]. For the evaluation of one-loop integrals OpenLoops employs the Collier library [33][34][35][36] and, alternatively, CutTools [37,38] together with the OneLOop library [39]. While we do not apply the resonance-aware method [22], the Powheg-Box-Res framework allows us to make use of new technical features, such as the automated implementation of scale variations , a Rivet [40] interface and the option to unweight events partially. This ttbb generator will be soon publicly available on the Powheg-Box webpage [23].

Powheg methodology
In the following, we briefly review the Powheg method [6,7] with emphasis on the separation of radiation into singular and finite parts. In this context we discuss technical subtleties that arise in the case of multi-scale processes like pp → ttbb.
The master formula for the description of NLO radiation in the Powheg approach consists of two contributions, which arise from the splitting of real emission into singular (s) and finite (f) parts, Here R(Φ R ) should be understood as squared real-emission matrix element, and Φ R as the corresponding phase space. Similarly, Born and virtual contributions in the Born phase space are denoted as B(Φ B ) and V (Φ B ). The splitting (3.2) is implemented as is a damping function that fulfils F → 1 and F → 0, respectively, in the infrared and hard regions of phase space (see below).
The singular part of real radiation is resummed according to the Powheg formula where real emission is further split into FKS sectors [41], which isolate collinear singularities arising from individual emitters. In each sector, the emission phase space Φ R is factorised into the Born phase space Φ B and a one-particle radiation phase space Φ rad through an appropriate FKS mapping, The term within squared brackets in (3.4) generates the hardest radiation according to an emission probability R/B. The parameter k T,α = k T,α (Φ rad ) stands for the hardness of the radiated parton, and radiation harder than k T,α is excluded by means of corresponding Sudakov form factors, The term ∆(q cut ) in (3.4) represents the no-emission probability above the infrared cutoff q cut , and Sudakov form factors account for unresolved multiple emissions in a way that cancels infrared singularities while preserving the differential NLO cross sectionB(Φ B ) in the Born phase space. The latter is defined by integrating out the singular part of real radiation,B Here infrared cancellations between V and R s are controlled via FKS subtraction. The remaining finite part of NLO radiation is treated as in fixed-order calculations, Note that flux and symmetry factors as well as the convolution with PDFs are implicitly understood in (3.4) and (3.9).
Let us now come back to the details of the separation of the singular and finite parts of real emission in (3.3). Technically, the damping function F is implemented based on the kinematics of the actual FKS sector, i.e. (3.10) The default functional form of F in Powheg-Box [42,43] is is the usual factor that smoothly shifts the weight of real radiation from R s to R f when the hardness of the emission, k T,α , becomes of the order of h damp or higher. Note that the freely adjustable h damp parameter in Powheg plays an analogous role as the resummation scale µ Q in the MC@NLO method. This is because both parameters act as a k T -threshold that separates the radiative phase space into a hard region, which is described by fixedorder matrix elements, from a singular region, where large logarithms of soft and collinear origin are resummed to all orders by means of Sudakov form factors. More precisely, in the MC@NLO approach the factor R s,α /B and the terms ∆ in 3) is entirely determined by the matrix element, which also dictates the scale at which the shower starts emitting further partons. Thus the Powheg method has the advantage of being essentially independent on the shower starting scale. More generally, thanks to the fact that the first emission is completely independent of the parton shower, Powheg predictions are characterised by a rather mild sensitivity to systematic uncertainties associated with the parton shower.
In addition to the well-known h damp -dependent damping mechanism (3.12), the Powheg-Box also implements a theta function 5 of the form [42,44] where R α corresponds to the infrared (soft and collinear) approximation of the full matrix element. Schematically it has the factorised form with an FKS kernel K α (Φ rad ) and an underlying Born contribution B(Φ B ), whose kinematics is determined by the inverse of the mapping (3.6) in the actual sector α. By default, 4 Note that in the MC@NLO implementation of [12] the scale µQ is also taken as starting scale for the showering of so-called hard events, which represent the MC@NLO counterpart of (3.9). Instead, in Powheg such events are showered starting from the actual kT of the first emission. the cut-off parameter h bzd in (3.13) is set equal to 5. In this way, in the vicinity of IR singularities, where R α /R α → 1, radiative contributions are attributed to R s and resummed according to (3.4). On the contrary, when the real emission matrix element largely exceeds the IR approximation (3.14), the resummation of the full R/B kernel according to (3.4) is not well justified, and corresponding events are attributed to the finite remnant (3.9) through the theta function (3.13). In the standard Powheg-Box, and in Ref. [13], the damping function (3.11) is applied only to initial-state radiation. However, in the present ttbb generator we have extended it to all (massless or massive) final-state emitters, that have a FKS sector associated with it, in order to ensure a consistent resummation of QCD radiation off b-quarks.
was originally introduced in order to avoid possible divergences of R(Φ α )/B(Φ B ) due to so-called Born zeros, i.e. phase space regions where B(Φ B ) → 0. Such divergences cancel in theB/B ratio, i.e. they are not physical, and are not related to IR radiation. Corresponding Φ B regions should thus be attributed to the finite remnant. Otherwise they could lead to dramatic inefficiencies in the event generation. More generally, the damping factor (3.13) can play an important role also in case of multi-scale processes where the Born cross section involves enhancement mechanisms at scales well below the hard energy of the full process. Such enhancements can compete with the ones due to soft and collinear QCD radiation in a way that is somewhat analogous to Born zeros.
In the case of pp → ttbb, such effects can arise from the interplay of soft and collinear enhancements due to NLO light-jet radiation and to the generation of the bb system in regions with m bb m ttbb and/or p T,bb m ttbb . For example, let us consider a gg → ttbbg event with a gluon emission of ISR type. Its kinematics is generated starting from a gg → ttbb Born event through a mapping of type (3.6), which creates the required gluon recoil by boosting the final state of the gg → ttbb Born event in the transverse direction. The relevant boost factor, γ = 1/(1 − β 2 ) 1/2 , is determined by p T,j = p T,ttbb = γβE ttbb , where E ttbb is the ttbb energy of the bb system in the Born event. If we assume, for simplicity, that he gluon is emitted in the same azimuthal direction as the bb system in the Born event, then the bb transverse momentum of the radiative event becomes p T,bb = γ(p T,bb − βE bb ), where p T,bb and E bb are the bb transverse momentum and energy in the Born event. Thus, the FKS mapping can lead to a very significant reduction of p T,bb . More precisely, for radiative events with the effect of the FKS boost on the bb system amounts to Thus, since the bulk of the ttbb cross section is characterised by E ttbb E bb and γ ∼ 1, in the case of hard QCD radiation with˜ 1 the FKS mapping can lead to a drastic reduction of the p T of the bb system. As a result, in the region of small m bb , the ISR boost can enhance the R α (Φ α )/R α (Φ α ) ratio by up to a factor 6 (p T,bb /p T,bb ) 2 ∼ 1/˜ 2 . This violates the main assumption that justifies the Powheg formula (3.4) , which requires a sufficiently hard ttbb process as compared to the k T of NLO radiation. In particular, due to the sensitivity of the Born amplitude to scales of the order p T,bb ∼ (E bb /E ttbb ) p T,j p T,j , the factorisation formula (3.14) is not fulfilled.
Fortunately, this problematic behaviour emerges only in relatively hard regions of the Φ rad phase space. 7 Thus, as a remedy it is natural to shift such events into the finite remnant by means of the damping factor (3.13). In fact, in the case of ttbb production we have found that the h bzd -dependent cut plays an important role for the efficiency generation of Les Houches events (LHEs) as well as for a consistent scale dependence. Moreover, applying a large h bzd cut we have observed a significant enhancement of the QCD scale dependence. This can be attributed to the fact that scale variations in the soft term (3.4) are restricted to theB factor (3.8), where the unphysical distortions of the bb kinematics induced by the FKS mappings can jeopardise the natural cancellation of virtual and real contributions associated with a given Born configuration.
As discussed in Section 4, Powheg predictions for ttbb observables are rather stable with respect to variations of h bzd . Thus, in order to avoid an unphysical enhancement of the scale dependence, we have reduced h bzd from its default value of 5 to h bzd = 2. This guarantees a more reasonable consistency with the fixed-order scale dependence without shifting an excessive fraction of the cross section from dσ s to dσ f .

Input parameters, PDFs and scale choices
The predictions in Sections 4-5 are based on the following input parameters, scale choices and PDFs. Heavy-quark mass effects are included throughout using All other quarks are treated as massless in the perturbative part of the calculations. Since we use massive b-quarks, for the PDF evolution and the running of α S we adopt the 4F scheme. Thus, for consistency, we renormalise α S in the decoupling scheme, where top-and bottom-quark loops are subtracted at zero momentum transfer. In this way, heavy-quark loop contributions to the evolution of the strong coupling are effectively described at first order in α S through the virtual corrections. For the calculation of hard cross sections at LO and NLO, as well as for the generation of the first Powheg emission, we use the NNPDF30 nlo as 0118 nf 4 parton distributions [45] as implemented in the LHAPDFs [46] and the corresponding α (4F) S . 8 To assess 6 The 1/p T,bb dependence arises from the collinear singularity associated with the emission of the parent gluon in the topologies of type (d) in Fig. 5. 7 This holds also for similar issues due to final-state radiation. 8 More precisely, α PDF uncertainties we re-evaluate the weights of LHEs with 100 different PDF replicas, while using the nominal PDF set for parton showering.
Since it scales with α 4 S , the ttbb cross section is highly sensitive to the choice of the renormalisation scale µ R , and this choice plays a critical role for the stability of perturbative predictions. Following [8,11], we adopt a scale choice of the form with the scale-variation factor ξ R ∈ [0. 5,2]. This dynamic scale choice accounts for the fact that ttbb production is characterised by two widely separated scales, which are related to the tt and bb systems and are chosen as the geometric average of the respective transverse energies, The transverse energies E T,i = m 2 i + p 2 T,i are defined in terms of the rest masses m i and the transverse momenta p T,i of the bare heavy quarks. The scales (3.19) are computed according to physical kinematics, i.e. without projecting real emission events to the underlying Born phase space. The choice (3.18) is applied to all (N)LO matrix elements apart from the α S factor that results from the R/B ratio in (3.4). In that case α S is evaluated at the transverse momentum of the hardest Powheg emission, and that α S (k T,α ) factor is not subject to scale variations.
For the factorisation scale µ F we use 9 where ξ F ∈ [0.5, 2], and the total transverse energy of the ttbb system, H T , is computed in terms of bare-quark transverse momenta including also QCD radiation at NLO. Our nominal predictions correspond to ξ R = ξ F = 1, and to quantify scale uncertainties we take the envelope of the seven-point variation (ξ R , ξ F ) = (0.5, 0.5), (0.5, 1), (1, 0.5), (1, 1), (1, 2), (2, 1), (2,2). For the Powheg-Box parameters h bzd and h damp , which control the resummation of NLO radiation according to (3.11)-(3.13) as discussed in Section 3.1, we set h bzd = 2 (3.21) and Here the various E T,i are defined in the underlying Born phase space. To account for the uncertainties associated with these choice we apply the independent variations h bzd = 2, 5, 10 and h damp = H T /4, H T /2, H T , 1.5 m t , varying both parameters one at a time. 10 The above choices for µ R , µ F and h damp , as well as the employed PDFs correspond to the setup recommended in [11]. 9 This choice does not coincide with the scale µF = 1 2 i=t,t ET,i adopted in [8]. However, this difference has a rather minor impact on our predictions. 10 The choice h damp = 1.5 mt corresponds to the default setting used for inclusive tt production in ATLAS.

Parton shower settings and variations
By default, LHEs are showered with Pythia 8.2 using the A14 tune 11 , where ISR and FSR parameters as well as the MPI activity have been tuned in a single step using most of the available tt ATLAS data from Run 1 [47]. In the A14 tune, m b = 4.75 GeV and α (5F) S (M Z ) = 0.127, both for ISR and FSR, while in the default Monash tune α (5F) S (M Z ) = 0.13650. Since the shower evolution is implemented in the 5F scheme we shower events using 5F PDFs. Specifically, we choose the NNPDF30 nlo as 0118 5F PDFs. 12 The interplay between Powheg and Pythia is controlled by the scalup parameter, which describes the hardness of radiation in LHEs and may be taken as starting scale for Pythia. However, in order to avoid inconsistencies due to the fact that the Pythia evolution variable does not coincide with the definition of hardness in Powheg, we apply the following two-step procedure based on the PowhegHooks class. Instead of starting below scalup, we instruct Pythia to generate radiation up to the kinematic limit by setting pythia.readString("SpaceShower:pTmaxMatch = 2"); pythia.readString("TimeShower:pTmaxMatch = 2"); Then, to guarantee the correct ordering of emissions in Powheg and Pythia, we apply a veto on each Pythia emission that is harder than scalup according to the Powheg-Box definition of hardness. This is achieved by setting pythia.readString("POWHEG:nFinal = 4"); pythia.readString("POWHEG:veto = 1"); The remaining PowhegHooks settings are left to their default values.
At LOPS level we set the shower starting scale equal to H T /2 and vary it up and down by a factor two in order to assess the related uncertainty. At NLOPS, the shower starting scale is dictated by the kinematics of real emission matrix elements in the Powheg method. Thus, at variance with NLOPS predictions based on the MC@NLO method, Powheg predictions are free from uncertainties related to the choice of the shower starting scale.
In order to assess uncertainties due to the parton-shower modelling of g → bb splittings we vary the parameter TimeShower:weightGluonToQuark, which permits to select various optional forms of the g → QQ splitting kernel in Pythia 8. The default is option 4, which corresponds to the splitting probability [48] The factor (1 − δ) 3 , which suppresses the production of high-mass bb pairs, is derived from the H → gbb matrix 11 More precisely we use the Pythia 8.219 version and the specific A14 tune for NNPDFs, which is based on the NNPDF2.3 LO PDFs. 12 We have checked that showering LHEs with the NNPDF30 nlo as 0118 nf 4 (used for the NLO hard cross section) or the NNPDF30 lo as 0130 nf 4 instead of the 5F NNPDFs does not induce any significant change in our results. element by interpreting m H as the mass of a gluon dipole, m dipole . Omitting the factor (1 − δ) 3 in (3.23) corresponds to option 2 and results in a DGLAP splitting probability of type γ * → bb with mass effects. More precisely, in option 2, g → bb splittings are generated based on massless kinematics, and the r b mass correction is implemented through reweighting, while massive kinematics is restored through momentum reshuffling. Option 3, which implements massive DGLAP splittings in a more realistic way, involves an additional (1+δ)/(1−δ) factor that leads to a significant enhancement of the g → bb rate. This option is excluded by LEP/SLC data and also by direct measurements of tt+b-jet production [49]. Finally, option 1 corresponds to option 2 with r b = 0 and yields very similar results. Thus, for the assessment of g → bb shower uncertainties we will compare options 4 and 2.
In addition to the functional form of the heavy-quark splitting kernel we also vary the scale of α S in the parton shower. To this end, we set TimeShower:weightGluonToQuark to 6 and 8, which corresponds to options 2 and 4 with α S (p 2 T ) replaced by α S (m 2 bb ) in the heavy-quark splitting kernel (3.23). Moreover, using TimeShower:renormMultFac, we vary α S (p 2 T ) → α S (ξp 2 T ) with prefactors ξ = 0.1, 1, 10 both for options 2 and 4. This latter variation is applied to all final-state QCD splittings, i.e. also splittings of type g → gg, q → qg, etc.

Comparisons against alternative generators
In order to assess systematic uncertainties related to the parton shower and the matching scheme, in Sections 4.2-4.3 we compare Powheg+Pythia predictions of ttbb production against corresponding predictions generated with Powheg+Herwig and with Sherpa [8]. The Powheg+Pythia and Sherpa generators of ttbb production are also compared against corresponding generators of inclusive tt production in the 5F scheme 13 .
In the case of Herwig In the case of Sherpa we use version 2.2.4 with its default tune 14 for Sherpa's dipole shower [52]. The relevant one-loop matrix elements are computed with OpenLoops, and 13 In the case of Powheg, the well known hvq generator [25] is used. 14 More precisely, in order to be consistent with the Sherpa 2.1 benchmarks presented in Ref. [11] we have used the shower recoil scheme proposed in Ref. [51], which was the default in Sherpa 2.1. This corresponds to setting CSS KIN SCHEME=0, while CSS KIN SCHEME=1, which became the new default in Sherpa 2.2, leads to slightly more significant differences with respect to the tt + b-jet predictions of Ref. [11]. More precisely, comparing Powheg+Pythia against Sherpa 2.2 with the new recoil scheme we observe differences at the level of 10% in the ttbb cross section and up to about 40% in the light-jet pT spectrum. Such differences are well consistent with QCD scale variations. For comparison they are three times smaller with respect to the differences between Sherpa+OpenLoops and Madgraph5aMC@NLO in Ref. [11]. Note that CSS KIN SCHEME acts only on the second and subsequent shower emissions, i.e. it does not affect the SMC@NLO matching procedure. matching to the parton shower is based on the Sherpa implementation [9] of the MC@NLO method [10], dubbed SMC@NLO. As for the hard cross section we use the same input parameters, PDFs and scale settings as specified in Section 3.2 for the case of Powheg. Moreover, as motivated in Section 3.1 we identify the resummation scale µ Q in Sherpa with the h damp parameter in Powheg, i.e. we set µ Q = H T /2. In the Sherpa simulation the NNPDF30 nlo as 0118 nf 4 PDF set is used throughout, i.e. also for paton showering.
For Powheg and Sherpa simulations of inclusive tt production we use the same setup as for the corresponding ttbb generators, with the only exceptions being the QCD scales, µ R = µ F = 0.5 E T,t E T,t , and the choice of the NNPDF30 nlo as 0118 nf 5 PDF set. In this setup the inclusive NLO cross section amounts to σ tt = 815 pb, which is only 2% below the NNLO prediction of 832 +45 −50 pb [53].

Simulations with stable or decayed top quarks
In Sections 4-5 predictions for ttbb production are presented both for the case of stable top quarks and with spin-correlated top decays. Simulations with stable top quarks permit to avoid the combinatorial complexity that results from the presence of four b-quarks in decayed ttbb events. In this way one can focus on the production of the bb pair that is governed by QCD dynamics and which represents the main source of theoretical uncertainty in pp → ttbb. Moreover, results with stable top quarks can be compared to the benchmarks of Refs. [8,11]. Since top quarks do not hadronise, when we switch off top decays we disable hadronisation and, following Ref. [8,11], we also deactivate multi-parton interactions (MPIs) and QED radiation in the parton shower. This is achieved by setting pythia.readString("6:mayDecay = off"); pythia.readString("-6:mayDecay = off"); pythia.readString("SpaceShower:QEDshowerByQ = off"); pythia.readString("SpaceShower:QEDshowerByL = off"); pythia.readString("TimeShower:QEDshowerByQ = off"); pythia.readString("TimeShower:QEDshowerByL = off"); pythia.readString("PartonLevel:MPI = off") pythia.readString("HadronLevel:All = off"); For the case of decaying top quarks we show results both with hadronisation and MPI switched off or on, while the QED shower is always activated and hadrons are kept stable throughout. For the implementation of spin-correlated decays in the Powheg framework we follow the approach of Ref. [54], which has already been employed in the Powheg-Box framework in Refs. [25,55,56]. More precisely, we use resonant tree matrix elements for the full 2 → 8 Born processes qq/gg → t(→ bij)t(→bkl)bb, where ij and kl stand for the leptons or quarks from W decays, and corresponding 2 → 9 processes with an additional external gluon at the level of the pp → ttbb sub-process. In the 2 → 8(9) matrix elements we include only topologies with two intermediate top resonances. This accounts for spin correlations as well as for off-shell effects associated with the top and the W propagators. Technically, top decays are generated starting from on-shell ttbb events with a veto algorithm based on the ratio between 2 → 8(9) matrix elements and corresponding 2 → 4(5) matrix elements for the underlying pp → ttbb(+jet) process.
As additional input parameters for top decays we use [57] M W = 80.385 GeV, the total widths Γ t = 1.329 GeV, Γ W = 2.089 GeV , (3.25) and the branching ratios where we assume a 100% branching ratio for t → bW decays. For the total W -boson branching ratios into leptons and hadrons we use the values [57] BR W →had = 0.675, which include state-of-the-art higher-order corrections.

Jet observables and acceptance cuts
For the reconstruction of jets we use the anti-k T [58] algorithm with R = 0.4. We select jets that fulfil p T > 25 GeV, |η| < 2.5, (3.29) both for the case of light jets and b-jets. At parton level, we define as b-jet a jet that contains at least a b-quark, i.e. jets that contain a bb pair arising from a collinear g → bb splitting are also tagged as b-jets. At particle level, i.e. when hadronisation is switched on, we tag as b-jets those jets that are matched to a B-hadron using the ghost method as implemented in FastJet [59]. When studying ttbb production with stable top quarks, in Sections 2 and 4, we categorise events according to the number N b of b-jets that do not arise from top decays and fulfil the acceptance cuts (3.29). For the analysis of cross sections and distributions we consider an inclusive selection with N b ≥ 1 and a more exclusive one with N b ≥ 2. We refer to them as ttb and ttbb selections, respectively.
In Section 5 we present predictions for ttbb production with top-quark decays in the dilepton channel. In this case we require two oppositely charged leptons, = e or µ, with p T, > 20 GeV, |η | < 2.5. Charged leptons are dressed with collinear photon radiation within a cone of radius 0.1. We do not apply any cut on missing transverse energy. Jets are defined as for the case of stable top quarks, and we select events with at least four b-jets that fulfill the acceptance cuts (3.29). TeV and their ratios in the phase space regions with N b ≥ 1 (ttb) and N b ≥ 2 (ttbb) b-jets as well in the region m b1b2 > 100 GeV of the ttbb phase space (ttbb 100 ). Nominal fixed-order predictions at LO and NLO accuracy are compared to corresponding LOPS and NLOPS predictions of Powheg+Pythia. Also cross sections at LHE level are reported. Uncertainties correspond to the envelope of the 7-point factor-two variations of µ R and µ F .

Predictions for ttbb production with stable top quarks
In this section we present numerical predictions for pp → ttbb at √ s = 13 TeV in the 4F scheme. The presented results have been obtained with Powheg+OpenLoops using the setup of Section 3. Top quarks are kept stable throughout as specified in Section 3.5, and we study cross sections and distributions in the inclusive ttb phase space with N b ≥ 1 b-jets, as well as in the ttbb phase space with N b ≥ 2.

NLOPS predictions with perturbative uncertainties
In this section we compare (N)LO and (N)LOPS predictions focusing on NLO and matching effects as well as perturbative and PDF uncertainties. Table 1 presents cross sections in the ttb and ttbb phase space, as well as in the presence of an additional cut, m b 1 b 2 > 100 GeV, on the invariant mass of the two hardest b-jets. At fixed order we find perfect agreement with the NLO results of Ref. [11]. The various phase space regions feature similar NLO uncertainties, around 25-30%, while corresponding LO scale variations are roughly a factor two larger. Both at LO and NLO, scale uncertainties are strongly dominated by µ R variations. The large σ ttb /σ ttbb ratio, which exceeds a factor 5, reflects the appearance of large logarithms of m b when a b-quark becomes unresolved. As shown in Section 2.3, such logarithms are mainly due to FS g → bb splittings. Thus the use of 4F PDFs, where ln(m b /Q) effects of IS origin are not resummed in the PDF evolution is well justified. Note also that ln(m b /Q) effects in σ ttb are present already at LO. Thus they do not jeopardise the convergence of the perturbative expansion. In fact, σ ttb /σ ttbb turns out to be very stable with respect to NLO corrections. The same hold for σ ttbb /σ ttbb 100 .
At variance with [8], where LO calculations were performed using LO PDFS and the corresponding value of α S , here, in order to obtain a more realistic picture of the convergence of the α S -expansion, we use NLO inputs throughout. 15 This approach increases the NLO K-factors from 1.15-1.25 [8] to 1.80-1.95. This observation raises some concerns regarding the possible presence of significant higher-order corrections beyond NLO and calls for a better understanding of the origin of the large K-factor at NLO. This question as well as the search for possible improvements is deferred to future studies.
Comparing fixed-order (N)LO cross sections against (N)LOPS ones we find that matching and showering effects are almost negligible in σ ttb , while in the case of σ ttbb they slightly exceed 10%, and in the Higgs-signal region, m b 1 b 2 > 100 GeV, they approach 30% . As pointed out in Ref. [8], such effects can be understood in terms of tt + 2b-jet production via double g → bb splittings. In practice, one of the b-jets results from a g → bb splitting in the ttbb matrix element, while the second one is created by the parton shower via a further g → bb collinear splitting. This interpretation is confirmed by the fact that the enhancement at hand is not present in the LHE-level cross sections presented in Table 1. In fact, double splittings are generated only at NLOPS level through parton showering. Double-splitting enhancements in Table 1 behave in a qualitatively similar way as in Refs. [8,11], but their size turns out to depend on the employed NLOPS generator. As compared to Ref. [11], we observe that the NLOPS/NLO correction to σ ttbb in Table 1 (+12%) is twice as large as in Sherpa (+6%), very close to Powhel 16 (+13%) and well below the prediction of For what concerns scale variations, in Table 1 we see that their impact at NLOPS tends to be 5-10% higher as compared to fixed-order NLO. This is consistent with the behaviour of Madgraph5aMC@NLO and Powhel in Ref. [11], while Sherpa features a significantly lower scale uncertainty. Such differences may be an artefact of the incomplete implementation of scale variations in the various NLOPS tools. In the case of Powheg, as anticipated in Section 3.1 we have found that increasing h bzd can lead to unphysical enhancements of the scale uncertainty. This effect is mostly visible in the ttbb phase space, where the maximum scale variation amounts to +40% for h bzd = 2 and grows up to +45% and +54% when setting h bzd = 5 and 50, respectively. Based on these observations, as default for our ttbb simulations we have set h bzd = 2. This choice guarantees a decent consistency with fixed-order scale variations without altering the matching procedure in a drastic way. In particular, when h bzd is reduced from its standard Powheg-Box value of 5 down to 2, we have checked that the fraction of the tt + b-jet cross section that is shifted from the singular part (3.4) to the finite remnant (3.9) amounts to only 10-20%. This holds for all considered distributions in the ttb and ttbb phase space.
LO inputs for the LO cross section results is a very strong dependence on the LO value of αS. The latter can depend very strongly on the employed PDF set. In particular, the two existing NNPDF 4F LO sets, which correspond to αS(MZ ) = 0.118 and αS(MZ ) = 0.130, can result in factor 1.5 ambiguity in the ttbb K-factor. We also note that the local K-factor in the Powheg matching formula, i.e. theB/B ratio in (3.4), is computed using NLO input throughout. 16 We note in passing that for other observables, especially in the ttb phase space, results in this paper can deviate more significantly from the Powhel predictions of Ref. [11]. This can be attributed to differences in the Pythia settings and, most importantly, to the fact that the Powhel generator used in Ref. [11] was based on 5F ttbb matrix elements with m b = 0 and made use of technical generation cuts in order to avoid collinear singularities from g → bb splittings. This limitation has now been overcome by upgrading the Powhel ttbb generator to the 4F scheme [13]. that result from the interplay of real-emission matrix elements and g → bb parton-shower splittings. Thus they feature an enhanced scale dependence.
For kinematic distributions that are inclusive with respect to NLO QCD radiation, NLOPS scale variations have a minor impact on shapes and amount essentially to a normalisation shift, similar to what observed at the level of the ttb and ttbb cross sections. In contrast, in the case of the light-jet p T spectra, scale variations increase from about 30% in the soft region up to 100% in the hard tails. This is consistent with the fact that such observables are only LOPS accurate and depend on α 5 S (µ R ). The effect of PDF variations is clearly subleading as compared to scale uncertainties and has little impact on shapes.
Comparing (N)LOPS predictions to the respective fixed-order (N)LO results, we observe that matching and shower effects remain almost negligible also at the level of distributions in the ttb phase space. As for the ttbb region, the NLOPS effects of order 10% observed in σ ttbb turn out to be quite sensitive to the kinematics of b-jets. In particular, as expected from the QCD dynamics of double g → bb splittings [8], the most pronounced effects are observed in the tails of the m b 1 b 2 and ∆R b 1 b 2 distributions, where the NLOPS/NLO ratio approaches a factor two. In the Higgs signal region, m b 1 b 2 ∼ 125 GeV, the NLOPS enhancement is around 1.25 and well consistent with Ref. [8].
Comparing fixed-order NLO predictions to LO ones we find that, in spite of the fairly large K-factors observed in Table 1, the shapes of distributions turn out to be quite stable with respect to higher-order QCD corrections. In the case of (N)LOPS predictions, the situation is different, especially for the shape of the light-jet p T spectra, which receives significant NLO distortions. This is not surprising, since at LOPS the light-jet p T is entirely generated by the parton shower. Thus the NLOPS/LOPS ratio should be regarded as a LO matrix-element correction to the parton-shower approximation, rather than a NLO correction in the perturbative sense.
Significant differences between NLOPS and LOPS shapes are observed also in the m b 1 b 2 and ∆R b 1 b 2 distributions. Since the respective NLO and LO shapes are very similar, this behaviour can be attributed to the parton shower. More precisely, it can be understood as a side effect of the above-mentioned NLOPS/LOPS correction to the light-jet p T spectra, which is converted into a double-splitting effect by g → bb splittings inside the light jet.

Shower uncertainties
In Figures 9-10 we study the sensitivity of (N)LOPS predictions to parton-shower and matching uncertainties for the same observables considered in Section 4.1.
The ratios displayed in the upper frames illustrate the net effect of parton showering by comparing full NLOPS predictions against results at LHE level. In addition, to assess parton-shower uncertainties, NLOPS predictions based on Pythia are compared to the corresponding results obtained with Herwig.
In the ttb phase space, apart from a mild distortion of the light-jet spectrum, the net effect of parton showering is essentially negligible. In contrast, in the ttbb phase space it increases the cross section by about 5% and tends to grow in the tails of distributions. The most sizable shower effects are observed in the m b 1 b 2 and ∆R b 1 b 2 distributions, where they reach up to 50-100%. This behaviour is well consistent with the enhancement of the NLOPS/NLO ratio observed in Fig. 8, and the fact that it is driven by the parton shower provides further support to its interpretation in terms of double g → bb splittings.  Predictions and uncertainties as in Fig. 9.
In spite of the important role of parton showering, it is reassuring to observe that the sensitivity of NLOPS predictions to the choice of parton shower is very small. In fact, the typical agreement between results based on Pythia and Herwig is at the level of a few percent both in the ttb and ttbb selections. Sizeable deviations at the level of 20% are observed only when requiring more than three b-jets. As discussed in Section 3.1, the very mild sensitivity of Powheg predictions to the choice of parton shower is due to the fact that the first emission is completely independent of the parton shower in the Powheg approach.
The ratios shown in the central frames of Figures 9-10 illustrate (N)LOPS uncertainties related to the modelling of g → bb splittings and variations of α S in Pythia (see Section 3.3). At LOPS also variations of the shower starting scale (scalup) are shown.
The fact that ttbb 4F matrix elements populate the whole bb phase space restricts the effect of g → bb shower splittings to events with four or more b-quarks. Thus, only the cross sections with N b ≥ 3, 4 b-jets suffer from sizable shower uncertainties. Vice versa, all considered observables with ttb or ttbb cuts turn out to be very stable, with typical shower uncertainties of a few percent at NLOPS. This holds also for the observables that are most sensitive to double splittings, i.e. m b 1 b 2 and ∆R b 1 b 2 , the only exception being the tail of the ∆R bb distribution, where double-splitting effects can reach 50% of the NLOPS cross section, while g → bb shower uncertainties can reach 15%.
Predictions at LOPS depend also on the choice of the shower starting scale. This uncertainty is especially sizeable in the case of the light-jet spectrum, where scalup acts as a cutoff. A sizeable scalup dependence is visible also in the LOPS predictions for the p T -distributions of b-jets, which indicates that such observables are rather sensitive to QCD radiation. Let us recall that the scalup dependence disappears completely in NLOPS simulations based on the Powheg approach.
Ratios plotted in the lower frames of Figures 9-10 show the dependence of NLOPS predictions with respect to the choice of the h damp and h bzd parameters, which control the separation of the first emission into events of soft and hard type in the Powheg-Box framework (see Section 3.3). The h damp band is obtained by varying h damp = H T /4, H T /2, H T , 1.5m t with the value of h bzd fixed to 2, while the h bzd band is obtained by varying h bzd = 2, 5, 10 with fixed h damp = H T /2. Observables that are inclusive with respect to light-jet radiation reveal a remarkably small dependence, typically of the order of a few percent, on the choice of h damp and h bzd . Non-negligible but moderate uncertainties are found only in the light-jet spectra, which are enhanced by up to 20% when h bzd is increased from 2 to 10. Investigating simultaneous variations of h damp and h bzd (not plotted) we have found that the size of the h damp variation band is fairly stable with respect to the value of h bzd within the considered range.

Comparisons against other ttbb and tt generators
In Figures 11-12 we compare tt+b-jet predictions based on Powheg+Pythia and Sherpa. This comparison is done both for (N)LOPS pp → ttbb generators in the 4F scheme and for corresponding generators of inclusive tt production in the 5F scheme. Specifically, in the case of Powheg we use hvq [25]. As detailed in Section 3.4, input parameters, QCD scales and matching parameters are chosen as coherently as possible across all generators. In this spirit, the parameter h damp = H T /2 in Powheg is identified with the resummation scale µ Q in the SMC@NLO framework of Sherpa. Instead, for what concerns the parton showers we simply use standard settings, i.e. we do not try to improve the agreement between generators by tuning the Pythia and Sherpa showers.
The ratios in the upper frames of Figures 11-12 show Powheg pp → ttbb predictions normalised to corresponding Sherpa predictions at LOPS and NLOPS accuracy. The bands describe the combination in quadrature of all matching and shower uncertainties 17 in Powheg+Pythia (referred to shower uncertainties in the following), while only nominal Sherpa predictions are considered in the ratios. Comparing LOPS predictions gives direct insights into the different modelling of radiation in Pythia and Sherpa. For observables that are inclusive with respect to jet radiation we find deviations between 10-40% and comparably large shower uncertainties. In contrast, in the jet-p T distributions the LOPS predictions of Pythia are far above the ones by Sherpa, with differences that can reach a factor 2.5 in the tails. These differences are perfectly consistent with LOPS shower uncertainties, which are dominated by variations of the Pythia starting scale.
Moving to NLOPS reduces the direct dependence on the parton shower. At the same time, differences between the Powheg and SMC@NLO matching methods come into play. In practice, at NLOPS we observe a drastic reduction of shower uncertainties, especially in the light-jet and b-jet p T -distributions. Also the differences between Powheg and Sherpa become very small at NLOPS. The ttb and ttbb cross sections agree at the percent level, and differential b-jet observables deviate by more than 5% only in the tails of the m b 1 b 2 and ∆R b 1 b 2 distributions. Even the light-jet spectra in the ttb and ttbb phase space deviate by less than 10-20% up to high p T , in spite of the limited formal accuracy (LOPS) of such observables. In the light of these results, NLOPS theoretical uncertainties related to the matching scheme and the parton shower seem to be well under control in pp → ttbb. In particular, their impact appears to be clearly subleading as compared to QCD scale uncertainties.
In the central frames of Figures 11-12 we compare (N)LOPS generators of inclusive tt production based on Powheg+Pythia and Sherpa. In this case, the g → bb final-state splittings that give rise to tt + b-jet signatures are entirely controlled by the parton shower. At LOPS, also the parent gluon that splits into bb is generated by the parton shower. Nevertheless, the ttb and ttbb LOPS cross sections predicted by Powheg and Sherpa deviate by less than 30%-40%. Instead, as expected, the shapes of tt + b-jet observables vary very strongly, and in all considered light-jet and b-jet distributions Pythia results exceed Sherpa ones by a factor of two and even more. This excess is well consistent with the estimated LOPS shower uncertainties. At NLOPS, only g → bb splittings are controlled by the parton shower, while the emission of their parent gluon is dictated by LO matrix elements. Consequently, we observe a drastic reduction of shower uncertainties as compared to LOPS. Also the differences between Powheg and Sherpa are largely reduced at NLO, nevertheless they remain quite significant in various distributions.  To provide a more complete picture of the uncertainties of inclusive tt simulations, in the lower frames of Figures 11-12 we compare Powheg+Pythia generators of inclusive tt production and ttbb production. Shower uncertainties are shown only for the tt generator. At LOPS, the tt generator is strongly sensitive to the modelling of pp → ttg through initialstate gluon radiation in Pythia. As a result, the tt generator overestimates the ttb and ttbb cross sections by about 90% and 50%, respectively. This excess is strongly sensitive to scalup, and in the p T -distributions it is confined to the regions below 100-200 GeV, while the tails are strongly suppressed. Also the m b 1 b 2 and ∆R b 1 b 2 distributions feature strong shape differences as compared to LOPS ttbb predictions.
Such differences go down significantly at NLOPS. The ttb and ttbb cross sections predicted by the tt generator overshoot ttbb results by only 15-20%, and also b-jet observables feature an improved agreement with ttbb predictions. Nevertheless, in b-jet observables we find quite significant shape differences, especially for the m b 1 b 2 and ∆R b 1 b 2 distributions, and shower uncertainties remain far above the ones of the ttbb generator (see upper frame). As for the light-jet spectra, tt predictions turn out to lie above ttbb ones by about a factor of two in the tails. In principle, with the help of parton shower tuning NLOPS tt generators may be amenable to a reasonable description of inclusive tt+b-jet observables. However, in the light of the above results it should be clear that NLOPS ttbb generators are mandatory in order to achieve an acceptable level of shower systematics.

ttbb production with top-quark decays
In this section we present NLOPS results of the Powheg+Pythia ttbb generator with leptonic top-quark decays. More precisely we consider final states with oppositely charged leptons and/or muons. By default hadronisation and MPI are deactivated in Pythia, and their effect is shown separately. As detailed in Section 3.5, our implementation of top decays is based on resonant pp → ttbb → 2 2νbb(+j) matrix elements, where spin correlations are consistently taken into account.
Top-quark decays are due to weak interactions and, up to small corrections of O(Γ t /m t ), their effect factorises with respect to ttbb production. Thus, while they strongly increase the complexity of tt + b-jet events, top decays are not expected to interfere with the QCD dynamics of pp → ttbb in a significant way. In order to verify this hypothesis, in Fig. 13 we compare NLOPS pp → ttbb simulations with stable and decayed top quarks. To this end, based on Monte Carlo truth, all reconstructed jets are split into two subsets associated with ttbb production and top decays. Specifically, jets that contain a parton originating from showered top-decay products are attributed to top decays, otherwise to ttbb production. 18 At the level of top decays we require two b-jets and two charged leptons within the acceptance cuts (3.29)-(3.30), while for the "reconstructed" ttbb system we consider the same cuts and observables as for the case of stable top quarks. In order to mimic the leptonic branching ratio and the efficiency of acceptance cuts on top-decay products, the normalisation of the ttbb simulation with stable top quarks is adapted to the predictions   Figure 13. Distributions in the b-jets of the reconstructed ttbb system for pp → tt + b-jets with dileptonic top decays at √ s=13 TeV. Inclusive number of additional b-jets (a), distribution with ttb cuts in the p T of the first b-jet (b) and distributions with ttbb cuts in the p T of the first b-jet (c) and light-jet (d) as well as in the invariant mass (e) and ∆R (f) of the first and second b-jet. All results are based on Powheg+Pythia with hadronisation and MPI switched off. The ratio corresponds to NLOPS predictions with Pythia decays (ttbb+PSdecay) or stable top quarks (ttbb) normalised to corresponding ones with spin-correlated decays (ttbb+decay). Top-decay products are subject to acceptance cuts, while predictions with stable top quarks are normalised to ttbb+decay ones at the level of the ttbb cross section.
with decayed top quarks. This is done at the level of the ttbb cross section through a constant normalisation factor.
As shown in Fig. 13, tt + b-jet observables with stable top quarks and reconstructed top decays turn out to agree quite well: b-jet cross sections and distributions deviate by only 5-10%, and also in the light-jet p T -distribution decay effects hardly exceed 10%. These differences can be understood as indirect effect of the acceptance cuts on top-decay products, which result from the correlation between the kinematics of the tt system and the additional jets.
Keeping in mind that realistic b-jet observables consist of a combinatorial superposition of b-jets from ttbb production and from top decays, the fact that Monte Carlo truth  Figure 14. Predictions for pp → tt + b-jets at √ s=13 TeV after leptonic top-quark decays. Four b-jets and two leptons within acceptance are required without any distinction between b-jets from ttbb production and decay. Distributions in the inclusive number of b-jets (a) and in the p T of the first (b), second (c) and third (d) b-jet. All results are based on Powheg+Pythia, and in the lower frame nominal NLOPS predictions with spin-correlated decays without (ttbb+decay) and with hadronisation (ttbb+decay+HAD) and multi-parton interactions (ttbb+decay+HAD+MPI) are compared to corresponding ones with Pythia decays (ttbb+PSdecay). acceptance cuts on top decays have only a minor effect on the production of b-jets suggests that the essential features observed in pp → ttbb production, such as double-splitting effects, are expected to show up also in the presence of top decays.
In order to assess the importance of spin correlations, in Fig. 13 we also compare spin-correlated top decays to isotropic decays generated by Pythia. At the level of reconstructed ttbb observables this comparison does not reveal any significant effect of spin correlations.
A more realistic analysis of ttbb production and decay is presented in Figures 14-15 without any distinction between ttbb production and decay.
Comparing spin-correlated and isotropic top decays, in b-jet observables we find no significant deviation, and significant spin-correlation effects show up only in the azimuthal correlation of the two charged leptons.
In Figures 14-15 we also assess the relative impact of hadronisation and multi-parton interactions (MPI). It turns out that b-jet observables are very stable with respect to hadronisation, with differences between parton and hadron level that do not exceed the few percent level. The same holds for MPI effects.
The above results indicate that insights on the QCD dynamics of ttbb production gained through studies with stable top quarks at parton level should hold true also in the presence of top decays and hadronisation.

Summary and conclusions
Searches for ttH production in the H → bb channel call for a precise theoretical description of the irreducible tt + b-jet background. To shed light on the QCD dynamics that governs this nontrivial multi-scale process, in the first part of this paper we have analysed the relative importance of the various mechanisms that lead to the radiation of b-quarks off pp → tt events. To this end we have compared the role of pp → ttbb topologies involving initial-state and final-state g → bb splittings. Using a naive diagrammatic splitting, as well as gauge-invariant collinear approximations, we have demonstrated that the tt + b-jet cross section is strongly dominated by b-jet production via final-state g → bb splittings. This holds both for phase space regions with two or only one resolved b-jets. These findings support the usage of NLOPS generators based on pp → ttbb matrix elements in the fourflavour scheme, while we have pointed out that tt + b-jet predictions based on tt+multi-jet merging rely very strongly on the parton-shower modelling of g → bb splittings.
Motivated by these observations we have introduced a new pp → ttbb Powheg generator in the 4F scheme. This tool is based on the Powheg-Box-Res framework, and all relevant matrix elements are computed with OpenLoops. When applied to a multi-scale process like pp → ttbb, the Powheg method can lead to subtle technical issues. In particular, we have pointed out that the FKS mappings which generate the recoil associated with the first Powheg emission can enhance the amplitude of the underlying ttbb Born process in a way that leads to anomalously large weights as compared to the behaviour expected from the factorisation of soft and collinear radiation. Fortunately, such anomalies arise only from events with finite transverse momenta and not in the soft and collinear limits. Moreover, the Powheg-Box framework disposes of a mechanism that automatically attributes such events to the so-called finite remnant, where QCD radiation is handled as in fixed-order NLO calculations. This mechanism, which is controlled by the h bzd parameter in (3.13), plays an important role for the efficiency of event generation. Moreover, it permits to avoid artefacts that can result from the application of QCD factorisation and resummation far away from their validity domain.
We have discussed predictions of the new Powheg generator and theoretical uncertainties for various tt+b-jet cross sections and distributions at the 13 TeV LHC. At variance with previous studies, in order to provide a better picture of the perturbative convergence, we have evaluated QCD corrections using the same α S value and the same PDFs at LO and NLO. The resulting NLO K-factors turn out to be close to two, even if the renormalisation scale is chosen in a way that is expected to absorb large logarithms associated with the running of α S . The question of the origin of such large higher-order effects and the search for possible remedies, such as improved scale choices, deserve to be addressed in future studies.
Scale uncertainties at fixed-order NLO amount to 25-30% and are dominated by renormalisation-scale variations. At NLOPS they tend to increase in a similar way as in Madgraph5aMC@NLO, while in Sherpa they tend to decrease [11]. However this behaviour may be an artefact of the incomplete implementation of scale variations in NLOPS generators.
Comparing predictions at NLO, LHE and LOPS level reveals significant shower effects at the level of 10% in the ttbb cross section and up to 30% or more for the invariant-mass and ∆R distributions of b-jet pairs. These effects can be attributed to double g → bb splittings [8] and are qualitatively and quantitatively consistent with the findings of Refs. [8,11]. For the p T -distribution of light-jet radiation, the predictions of the new Powheg generator are quite close to fixed-order NLO and also quite stable with respect to variations of the parameters h damp and h bzd , which separate real radiation into singular and finite parts. This good stability is guaranteed by the h bzd -dependent mechanism mentioned above.
To assess pure shower uncertainties we have compared Powheg samples generated with Pythia 8 and Herwig 7. In addition, we have considered systematic uncertainties due to the modelling of g → bb splittings and the choice of α S in Pythia. At NLOPS, all shower uncertainties turn out to be rather small and clearly subleading with respect to QCD scale variations. As a further independent estimate of matching and shower uncertainties we have compared NLOPS ttbb generators based on Powheg+Pythia and Sherpa finding remarkable agreement both for tt + b-jet cross sections and distributions. We have also shown that matching and shower uncertainties increase considerably if NLO corrections are not taken into account. The same holds for NLOPS generators of inclusive tt production as compared to ttbb generators.
Finally, we have presented predictions for pp → ttbb with spin-correlated top decays. In this context we have show that hadronisation and MPI effects are almost negligible. Thus, the key features of the QCD dynamics of ttbb production at parton level are expected to hold true also at particle level after top decays.
The new ttbb Powheg generator will be made publicly available in the near future, and its application to experimental analyses may lead to significant steps forward in the understanding of the QCD dynamics of tt + b-jet production and in the control of the theoretical uncertainties that plague ttH(bb) searches.