NLO QCD predictions for $t\bar{t}b\bar{b}$ production in association with a light jet at the LHC

Theoretical predictions for ttbb production are of crucial importance for ttH measurements in the H->bb channel at the LHC. To address the large uncertainties associated with the modelling of extra QCD radiation in ttbb events, in this paper we present a calculation of pp->ttbbj at NLO QCD. The behaviour of NLO corrections is analysed in a variety of observables, and to assess theoretical uncertainties we use factor-two rescalings as well as different dynamic scales. In this context, we propose a systematic alignment of dynamic scales that makes it possible to disentangle normalisation and shape uncertainties in a transparent way. Scale uncertainties at NLO are typically at the level of 20-30% in fiducial cross sections, and below 10% for the shapes of distributions. The kinematics of QCD radiation is investigated in detail, including the effects of its recoil on the objects of the ttbb system. In particular, we discuss various azimuthal correlations that allow one to charaterise the QCD recoil pattern in a precise and transparent way. In general, the calculation at hand provides a variety of precise benchmarks that can be used to validate the modelling of QCD radiation in ttbb generators. Moreover, as we will argue, pp->ttbbj at NLO entails information that can be used to gain insights into the perturbative convergence of the inclusive ttbb cross section beyond NLO. Based on this idea, we address the issue of the large NLO K-factor observed in the ttbb cross section, and we provide evidence that supports the reduction of this K-factor through a mild adjustment of the QCD scales that are conventionally used for this process. The presented 2->5 NLO calculations have been carried out using OpenLoops 2 in combination with Sherpa and Munich.


Introduction
The associated production of top-and bottom-quark pairs at hadron colliders is an especially interesting process. From the theoretical point of view, it offers rich opportunities to investigate the dynamics of QCD in the presence of multiple scattering particles and energy scales. In particular, higher-order calculations of pp → ttbb raise non-trivial questions related to the mass gap between m b and m t , the choice of QCD scales, and the convergence of the perturbative expansion. Further strong motivation for a deeper understanding of ttbb production comes from its critical role as irreducible background to ttH production with H → bb at the LHC [1][2][3]. In this context, the modelling of pp → ttbb represents the main source of uncertainty in ttH(bb) measurements. Thus, improving the theoretical description of the ttbb background is of great importance for the sensitivity of ttH(bb) analyses at the High-Luminosity LHC [4]. Precise theoretical calculations for ttbb production are relevant also for direct experimental studies of this process, and recent measurements of the ttbb cross section [5][6][7] tend to exceed theory predictions by 30-50%.
At leading order (LO) in QCD, the ttbb cross section is proportional to α 4 S and suffers from huge scale uncertainties. Next-to-leading order (NLO) QCD calculations [8][9][10] reduce scale uncertainties to 20-30%, but the level of precision and the size of the corrections depend in a critical way on the choice of the renormalisation scale µ R . In this respect, in order to avoid an excessively large NLO K-factor, it was found that the value of µ R should be chosen in the vicinity of the geometric average of the energy scales of the tt and bb systems [10].
Calculations of pp → ttbb based on the five-flavour (5F) scheme [8][9][10], where b-quarks are treated as massless partons, are applicable only to the phase space with two resolved b-jets, while including b-mass effects in the four-flavour (4F) scheme makes it possible to obtain NLO predictions in the full tt + b-jets phase space [11], including regions where one b-quark is unresolved. The choice of the 4F scheme as opposed to the 5F scheme is also supported by the fact that initial-state g → bb splittings play a marginal role in tt + b-jets production, while the vast majority of b-jets originate via initial-state gluon radiation with subsequent g → bb splittings [12].
In order to be applicable to ttH(bb) measurements, NLO calculations of pp → ttbb need to be matched to parton showers. Nowadays, this can be achieved within various Monte Carlo frameworks [11, 12, 14? -16], using different matching methods and parton showers. Some of these generators are in good mutual agreement, but the overall spread of Monte Carlo predictions suggests that ttbb modelling uncertainties may significantly exceed the level of QCD scale variations, thereby spoiling NLO accuracy [17]. In this context, the uncertainties related to the modelling of extra QCD radiation that accompanies ttbb production play a dominant role.
Motivated by these observations, in this paper we present a NLO QCD calculation of ttbb production in association with one additional jet at the LHC. 1 Bottom-mass effects are included throughout using the 4F scheme. For the calculation of the required 2 → 5 one-loop amplitudes, which involve up to 25'000 diagrams in a single partonic channel, we use the latest version of the OpenLoops program [20], where scattering amplitudes are computed with the new on-the-fly reduction method presented in [21]. For the calculation of hadronic cross sections, OpenLoops 2 is interfaced with Sherpa [22][23][24][25] and, alternatively, with Munich 2 .
We discuss NLO predictions for pp → ttbbj at 13 TeV with emphasis on the assessment of perturbative uncertainties. To this end, we study conventional scale variations as well as different dynamic scales, and we point out that the effects of these two kinds of scale uncertainties are largely correlated. Based on this observation, we propose the idea of aligning dynamic scales to a natural scale, which can be defined using the maxima of the NLO variation curves as a reference. This prescription makes it possible to disentangle the effects of factor-two variations and dynamic scale variations in a way that provides a more transparent picture of normalisation and shape uncertainties.
To characterise the behaviour of QCD radiation in ttbb events, we consider kinematic distributions in the hardest light jet as well as recoil effects on the various objects of the ttbb system. To this end, we introduce azimuthal angular correlations that provide a transparent and perturbatively stable picture of recoil effects. Our NLO predictions for these and various other observables can be used as precision benchmarks to validate the modelling of QCD radiation in ttbb generators.
Finally, we exploit the calculation at hand to address the issue of the large NLO Kfactor observed in the integrated ttbb cross section [12]. In this respect, we note that the NLO corrections to pp → ttbbj correspond to the same order in α S as the NNLO corrections to inclusive ttbb production, i.e. O(α 6 S ). Thus they entail (partial) information on the behaviour of σ ttbb beyond NLO. Based on this idea, we use the ttbbj cross section at NLO to identify an optimal scale choice for the process pp → ttbb. The results of this analysis support a slight adjustment of the conventional ttbb scale choice, which results in a reduction of the ttbb K-factor and is also expected to attenuate NLO matching uncertainties.
The paper is organised as follows. In Sections 2-3 we outline the main ingredients of pp → ttbbj at NLO, and we document the employed input parameters, scale choices and acceptance cuts. In Section 4 we study the integrated cross sections and their scale dependence, and we check the safeness of our predictions with respect to Sudakov logarithms beyond NLO. Moreover, we propose the idea of disentangling shape and normalisation uncertainties by means of an alignment prescription for dynamic scales. Differential observables and shape uncertainties are presented in Section 5, where we also discuss recoil effects. Finally, in Section 6 we use ttbbj NLO predictions to identify an improved scale choice for inclusive ttbb production. Our main findings are summarised in Section 7.
2 Ingredients of the calculation 2.1 ttbbj production in the 4F scheme We investigate NLO QCD corrections to hadronic ttbbj production in the 4F scheme, i.e. we treat not only top quarks, but also bottom quarks with a finite mass throughout. The nonvanishing bottom mass renders g → bb splittings finite, which allows us to investigate also observables with unresolved b-jets and to apply the experimentally favoured definition of b-jets as all hadronic jets that contain at least one bottom (anti-)quark at the parton level. In particular, jets resulting from the clustering of b andb partons are considered b-jets as well. Accordingly, only hadronic jets that are constituted from light quarks q = d, u, s, c and gluons are considered light jets. In the 4F scheme, since no bottom (anti-)quarks appear as proton constituents, no further bottom (anti-)quarks are generated at NLO QCD. Thus all b-jets are generated by Feynman diagrams that contain exactly one bb pair. Input parameters, renormalization scheme and parton-distribution functions (PDFs) are chosen according to the 4F scheme, as detailed in Section 3.1.
The independent partonic channels contributing to pp → ttbbj at NLO are summarised in Table 1 Figure 1. Selected Born diagrams in the gg → ttbbg channel.
qq → ttbbg with q = d, u, s, c, where the latter gives rise to six quark-anti-quark and gluon-(anti)quark channels via permutations of q,q, g. Fig. 1 illustrates sample diagrams for the gluon-gluon channel, which is by far the dominant channel, with a contribution of about 77% (qg: 21%, qq: 2%). The dominant gg → ttbbg topologies are those where the bb pair is emitted from a g → bb splitting and the final-state gluon results from an initial-state g → gg splitting, while the tt pair is produced in a t-channel configuration. However, the impact of other topologies becomes quite prominent in certain phase-space regions, like e.g. at high invariant mass or ∆R separation of the bb system. See also Fig. 3 for the dominant gg → ttbb topologies.
At NLO in QCD, as usually the process receives contributions both from virtual and real corrections, which are separately divergent. To mediate these divergences between the different phase spaces, we rely on the dipole-subtraction formalism [26] in its extension to massive QCD partons [27].
The virtual corrections are constituted from both diagrams with a closed quark loop and diagrams that are generated from the LO ones by exchanging a virtual gluon between Figure 2. Selected gg → ttbbg one-loop diagrams (first row) and gg → ttbbgg real-emission diagrams (second row). any of the external or internal legs. Since all involved partons interact under QCD, the number of loop diagrams is more than a factor of 50 larger than the number of Born diagrams in the respective channels (see Table 1). While the quark-loop diagrams contain up to pentagon functions, the gluon-exchange diagrams require up to heptagon functions. Some sample diagrams for the latter are shown in Fig. 2 (first row), again for the dominant gg channel only.
The real-correction channels are constructed from the LO ones by either emission of another gluon or by the splitting of a gluon into a light qq pair. Including crossings of light partons between initial and final states, the channels listed in Table 1 result. In Fig. 2 (second row) we depict sample diagrams for the dominant all-gluon channel.

Tools and validation
The calculations presented in this paper have been performed with the automated frameworks Sherpa+OpenLoops and Munich+OpenLoops. Each of them completes the full chain of operations -from process definition to collider observables -that enter NLO QCD simulations at parton level.
In both frameworks virtual amplitudes are provided by OpenLoops 2 [20], the latest version of the OpenLoops matrix-element generator. One of the of main novelties of OpenLoops 2, which is used for the first time in the calculation at hand, is the combination of the original open-loop algorithm [28] with the recently proposed on-the-fly reduction method [21]. In this approach, the construction of loop amplitudes and their reduction to scalar integrals are combined in a single numerical recursion, which makes it possible to generate one-loop amplitudes in a way that avoids high tensorial ranks at all stages of the calculations. This results in a significant speed-up for multi-leg processes. Specifically, for the process at hand, the excellent CPU performance of OpenLoops 1 is further improved by a factor of three. For the treatment of numerical instabilities, the on-the-fly reduction algorithm is equipped by an automated stability system that combines analytic expansions together with a novel hybrid-precision system. The latter detects residual instabilities based on the analytic structure of reduction identities and cures them by switching from double (dp) to quadruple (qp) precision. Thanks to the local and highly targeted usage of qp, the typical qp overhead wrt dp evaluation timings is reduced from two orders of magnitude to a few percent.
The only external ingredients required by OpenLoops 2 are the scalar integrals [29], which are provided by the Collier library [30,31] by default, or by the OneLOop library [32] for exceptional qp evaluations. All amplitudes have been thoroughly validated against OpenLoops 1 [28], where the reduction is carried out based on the Denner-Dittmaier techniques [33,34] available in Collier or, alternatively, using CutTools [35], which implements the OPP method [36], together with the OneLOop library [32] for scalar integrals. Additionally, matrix elements have been cross-checked against the completely independent generator Recola [37,38].
All remaining tasks, i.e. the bookkeeping of partonic subprocesses, phase-space integration, and the subtraction of QCD bremsstrahlung, are supported by the two independent and fully automated Monte Carlo generators, Munich and Sherpa.
In Sherpa, tree amplitudes are computed using Comix [24], a matrix-element generator based on the colour-dressed Berends-Giele recursive relations [39], while one-loop amplitudes are provided by OpenLoops. Infrared singularities are cancelled using the dipole subtraction method [26,27], as automated in Comix, with the exception of K-and P-operators that are taken from the implementation described in [25]. Comix is also used for the evaluation of all phase-space integrals. Analyses are performed with the help of Rivet [40], which involves the FastJet package [41,42] to cluster partons into jets.
The parton-level generator Munich has been applied to several multi-leg processes at NLO QCD and EW accuracy, and as a key ingredient of the Matrix framework [43] it has been intensively applied to boson and diboson production at NNLO QCD. Munich provides a very efficient multi-channel phase-space integration with several optimizations for higher-order applications. All tree-level and one-loop amplitudes are supplied by Open-Loops through a fully automated interface. The implementation of the massive dipole subtraction formalism used in the present calculation has been extensively tested in the context of off-shell top-pair production in the 4F scheme [44], and very recently in the NNLO QCD production of tt pairs [45,46]. The implementation of phase-space cuts at generation and analysis level, as well as the event selection including jet algorithms are realized directly in Munich, without relying on external tools. Also the calculation of arbitrary (multi-)differential observables and the setting of dynamic scales are handled internally. Thereby Munich provides an independent cross-check of basically all remaining steps of the working chain.
Both tools have been validated extensively against each other for a representative selection of the results presented in this paper. All cross sections binned in b-jet and light-jet multiplicities (see Tables 3 and 4) have been validated at a precision level of 0.3% throughout for all scale choices. Moreover, most of the differential distributions presented in Section 5 have been cross-checked at the NLO level. For all compared observables we find agreement on the level expected from the statistical uncertainties of the two independent calculations.

Technical aspects and setup
In this section we specify the input parameters, PDFs, scale choices and acceptance cuts used in the calculations presented in Sections 4-6.

Input parameters, PDFs and scale choices
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. We present predictions for pp → ttbbj at √ s = 13 TeV. At LO and NLO we use throughout the 4F NNPDF parton distributions [47] at NLO, and the corresponding strong coupling. 3 PDF uncertainties are expected to play a rather subleading role, similarly as for pp → ttbb [12]. Thus we will base our predictions on the nominal PDF set, restricting our assessment of theoretical uncertainties to perturbative scale variations. 4

Renormalisation and factorisation scales
Since it scales with α 5 S , the ttbbj 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. Along the lines of [11,12,17], we adopt a dynamic scale that accounts for the fact that ttbb production is characterised by two widely separated scales, which are related to the tt and bb systems. To this end we define where 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, without applying any jet 3 More precisely we use the NNPDF30 nlo as 0118 nf 4 parton distributions, as implemented in LHAPDF [48], where α (4F) S (MZ ) = 0.112, which corresponds to α (5F) S (MZ ) = 0.118. 4 Using 100 replicas of the PDF set at hand we have checked that PDF uncertainties are at the level of 10% for the integrated ttbbj cross section and grow slowly with the pT of the various final-state objects, reaching at most 20% in the regions where event rates are suppressed by two orders of magnitude. We refrain from reporting further details on PDF uncertainties since they are strongly correlated to the ones observed in inclusive ttbb production [12] and thus only marginally relevant for the theoretical questions addressed in this paper.  Figure 3. Generic leading-order gg → ttbb topologies with final-state (a) and initial-state (b) g → bb splittings. The bulk of the ttbb cross section is dominated by topologies of type (a) with rather collinear splittings, while initial-state collinear splittings become important in the region of large ∆R bb [12].
algorithm at NLO. Also m 2 bb is defined in terms of the bare four-momenta of the (anti-)b quarks. As default choice for the renormalisation scale we adopt the geometric average of the various transverse energies and momenta of the ttbbj system, where the rescaling factor ξ R is typically varied in the range [0. 5,2]. This choice represents the natural generalisation of the widely used scale [11,12,17] for ttbb production. 5 The additional light-jet p T that enters (3.3) is defined using an auxiliary 6 k T -jet algorithm with R = 0.4, which is applied only to massless partons, i.e. excluding top and bottom quarks from the recombination, and is free from any restriction in p T and rapidity. In order to assess shape uncertainties, we consider three alternative dynamic scales with different kinematic dependences. The first one is defined as where the bb system enters through its invariant mass and its total transverse energy, This choice is motivated by the fact that m bb and E T,bb correspond to the virtualities of the QCD branching processes that dominate ttbb production, namely initial-state g → gg splittings followed by a final-state g → bb splittings (see Fig. 3). 5 The choices (3.3)-(3.4) are motivated by the fact that, to lowest order in the strong coupling, α 5 S (µ ttbbj ) = α 2 S (µ tt ) α 2 S (µ bb ) αS(pT,j) and α 4 S (µ ttbb ) = α 2 S (µ tt ) α 2 S (µ bb ). In this way, the coupling factors associated with the production of the tt and bb systems, plus the additional light jet for ttbbj production, are effectively evaluated at the corresponding characteristic scales, µ tt , µ bb and pT,j, avoiding large logarithms associated with the evolution of αS. 6 For the definition of physical observables a conventional anti-kT algorithm is used (see below).  Table 2. Naming scheme for phase-space regions with different inclusive multiplicities of b-jets ) and light jets (N j ≥ N min j ) that pass the acceptance cuts (3.11).
As further alternatives we consider two other dynamic scales, and which are defined in terms of the transverse energies of the jets, 8) and the total transverse energy, Here E T,j = p T,j for massless partons, and the sums run over all final-state QCD partons, always including NLO radiation and excluding only top quarks in the case of H T,jets . The factorisation scale µ F represents the maximum transverse momentum for initialstate radiation that is resummed in the PDFs. Thus it is typically chosen of the order of the halved hard-scattering energy. Following [11,12] we use 7 where ξ F ∈ [0.5, 2]. 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).

Jet observables and acceptance cuts
For the reconstruction of jets we use the anti-k T [49] algorithm with R = 0.4. We select b-jets and light jets that fulfil the acceptance cuts To be precise, the choice (3.10) agrees with the one used in [12] but differs from the choice µF = 1 2 i=t,t ET,i made in [11]. However, this difference has a minor impact on our predictions.
We define as b-jet a jet that contains at least one b-quark, i.e. jets that contain a bb pair arising from a collinear g → bb splitting are also tagged as b-jets. 8 Top quarks are kept stable throughout. When studying ttbbj production, we categorise events according to the number N b of b-jets and the number N j of light jets that fulfil the acceptance cuts (3.11). We always consider inclusive phase-space regions with N b ≥ N min b and N j ≥ N min j , and we label them as indicated in Table 2. For the analysis of cross sections and distributions, we always require one additional jet, and we consider an inclusive ttbj selection (N min b = 1) and a more exclusive ttbbj selection (N min b = 2).

Integrated cross sections for pp → ttbbj at 13 TeV
In this section we present numerical predictions for pp → ttbbj at √ s = 13 TeV in the 4F scheme. The results have been obtained with Sherpa+OpenLoops and Mu-nich+OpenLoops, using the setup of Section 3. Top quarks are kept stable throughout, and we study cross sections and distributions in the inclusive ttbj and ttbbj phase-space regions as defined in Table 2, applying the acceptance cuts (3.11). Perturbative scale uncertainties are assessed by means of seven-point factor-two scale variations and by comparing the various dynamic scales defined in Section 3.2.

Renormalisation scale dependence
A first picture of the perturbative behaviour of the ttbbj cross section is displayed in Fig. 4, where LO and NLO predictions based on the nominal scale choice (3.3) are plotted as a function of the renormalisation scale µ R . For each value of µ R , the effect of factor-two scale variations is illustrated through two bands, which correspond to the variation of µ R alone and the full 7-point variation of µ R and µ F . The results demonstrate that µ F variations play only a marginal role, especially at NLO. Thus, in the following we will focus on the µ R dependence. In Fig. 5 we plot the LO and NLO ttbbj cross section as a function of µ R for all four dynamical scales defined in (3.3),(3.5)-(3.7). For each choice the renormalisation scale is varied around its nominal value by a factor ξ R ∈ [1 /16, 16], while the factorisation scale is kept fixed at µ F = H T /2. The behaviour of the LO curves in Fig. 5 reflects the α Sdependence of the LO cross section, σ LO ∝ α 5 S , and corresponds essentially to the running of α S to the fifth power. To discuss the qualitative behaviour of Fig. 5 in more detail, let us consider the effect of µ R → ξµ R rescalings at LO, Here , and for small variations δξ,  NLO (red) with scales µ R = ξ R µ ttbbj and µ F = H T /2 are plotted as a function of the renormalisation scale factor ξ R . The main frame presents absolute predictions and corresponding 7-point factortwo variations of µ R and µ F , which are shown as uncertainty bands. The relative impact of such variations at LO and NLO is displayed in the two ratio plots, which show also a second uncertainty band corresponding to pure factor-two variations of µ R at fixed This is consistent with the LO curves of Fig. 5, where we observe that around the nominal scales (ξ = 1), reducing µ R by a factor 2 augments the LO cross sections by a factor close to 2 and vice versa, which corresponds to 5a 0 (µ) ∼ 1.
At NLO, the one-loop α S -counterterm cancels the ξ-dependence at O(α S ln ξ), resulting in a significant reduction of scale variations. In the vicinity of the nominal scales, factortwo variations go down to 10-25%, depending on the type of scale and the direction of the variation. As usually, the various NLO curves feature a stable point, which is located between ξ R = 1/2 and 1/3. In the region below the maximum, the NLO curves start falling quite fast, and between ξ R =1/6 and 1/8 they lead to negative cross sections. To avoid such a pathologic perturbative behaviour, the normalisation factors in the definition of µ T,tot and µ T,jets have been chosen in such a way that factor-two variations of the nominal scales do not enter the region below the NLO maximum. Concerning the NLO correction factors, K = σ NLO /σ LO , at ξ R 1 we find K ∼ 1.5 while the K-factor approaches one in the vicinity of the NLO maxima of the respective curves.
A striking feature of Fig. 5 is that, in spite of the rather different kinematic dependence of the various dynamic scales, the observed LO and NLO scale variations and K-factors have a fairly similar shape. In order to gain more insights into the origin of this behaviour, in the following we focus on the α S -dependence of the LO cross section. For the differential and integrated cross sections let us define where µ dyn (Φ) is a certain dynamic scale, Φ stands for the fully-differential final-state phase space, and the convolution with PDFs as well as acceptance cuts are implicitly understood. For the integrated cross section with dynamic scale µ dyn we can write where the result is expressed in terms of the α S -free cross sectionσ LO and the coupling factor α 5 S (μ dyn ), which corresponds to the average of α 5 S (µ dyn (Φ)). The above identity is nothing but a definition of the "average" scaleμ dyn , which depends both on the functional form of µ dyn (Φ) and on the applied phase-space cuts. Let us now consider scale variations, The effect of µ dyn → ξ µ dyn on α S µ dyn (Φ) can be expressed as where the α S (ξμ dyn ) prefactor on the rhs corresponds to a trivial rescaling ofμ dyn , while the term between square brackets depends on all moments of the distribution in ln (µ dyn (Φ)), Such moments may influence the scale dependence in a non-trivial way. However, their actual impact on the integrated cross section turns out to be marginal. This is due to the fact that QCD cross sections are typically dominated by phase-space regions with welldefined energy scales in the vicinity of the thresholds for producing massive final states and passing acceptance cuts. As a consequence, the distribution in ln (µ dyn (Φ)) is confined in the vicinity of its average value, ln(μ dyn ), and its higher moments are rather strongly suppressed. This implies for all n ≥ 1. More precisely, let us assume 9 that X n = a n 0 (ξμ dyn ) ln(µ dyn ) − ln(μ dyn ) n 1 , (4.9) for n ≥ 2. This implies that the expectation value of the rhs of (4.6) is dominated by the n = 0 term. Thus, under the above assumptions, the scale dependence of the LO cross sections (4.5) can be approximated as i.e. by a naive rescaling of α 5 S (μ dyn ). We have verified that this property is fulfilled with percent-level accuracy by all LO curves of Fig. 5. This means that, at the level of the integrated cross section, the various scales (3.3),(3.5)-(3.7) are equivalent to each other. More precisely, the scale dependence of σ LO with a given dynamic scale µ dyn,k can be related to the one of a fixed scale µ 0 by means of a constant rescaling which results into Therefore, as far as the scale uncertainty of σ LO and its normalisation are concerned, comparing different types of dynamic scales has no significant added value wrt simple ξ Rrescalings. For this reason, we advocate the usage of "aligned" dynamic scalesμ dyn,k , as defined in (4.11). In this way, the various dynamic scales have the same average value, and the uncertainties related to this common average value are accounted for by standard ξ R -rescalings, while the comparison of different scale definitions allows one to highlight the genuine kinematic effects that are inherent in their dynamic nature. Comparing aligned dynamic scales yields no significant effect at the level of integrated cross sections, but provides key information on shape uncertainties, since the average scalesμ dyn,k are sensitive both to the probed phase-space regions and to the detailed kinematic dependence of µ dyn,k (Φ). Vice versa, ξ-rescalings can be used to assess uncertainties in the normalisation of σ LO , whereas their impact on shapes is typically quite limited.
At LO, the above-mentioned alignment approach misses a crucial ingredient, namely a good criterion for the choice of a reference scale µ 0 . For pp → ttbbj, due to the very strong scale dependence induced by α 5 S , the choice of a well-behaved central scale is of crucial importance. At the same time, the presence of multiple scales, distributed from m b to m tt and beyond, renders this choice non-trivial. At NLO, a natural way of addressing the scalechoice problem is to exploit the presence of a characteristic scale given by the maximum of the NLO scale-dependence curves, µ max . The maximum itself is not necessarily an optimal scale choice, since its position is not guaranteed to be stable wrt higher-order corrections. Moreover, the flatness of the scale dependence around µ R = µ max tends to underestimate scale uncertainties. A more reasonable and conservative option, that will be adopted in this paper, is to set the central scale at µ R 2µ max . In this way, the range of factor-two scale variations extends over [µ max , 4 µ max ], covering the maximum itself as well as a relatively broad region where σ NLO is monotonically decreasing.
As observed in Fig. 5, the position of µ max depends on the choice of the dynamic scale. However, for reasons similar to those discussed above at LO, also NLO scale variations and the position of their maxima can be aligned via rescalings. This is not entirely obvious and does not work as precisely as in the LO case. The main reason is that NLO cross sections consist of two kind of contributions: Born and virtual parts, which are distributed in a similar way as dσ LO , and real-emission parts that can be distributed in a significantly different way. Moreover, dynamic scales can feature a different sensitivity to the kinematics of hard jet radiation, leading to genuinely new scale-dependence effects at NLO. For these reasons, the LO scale-dependence model (4.3)-(4.12) should be refined by splitting σ NLO into two parts with independent average scales. Nevertheless, for the process at hand and the scale choices (3.3),(3.5)-(3.7), it turns out that a single overall rescaling can already yield a good level of NLO alignment. This is illustrated in Fig. 6, where the dynamic scales (3.5)-(3.7) have been rescaled in such a way that the positions of the NLO maxima match the maximum of σ NLO (µ ttbbj ), which is located at 0.45 µ ttbbj , i.e. µ ttbbj is rather close to 2µ max . This alignment is achieved is due to the fact that the latter had already been placed on purpose about two times above the maximum, but without tuning their position in a precise way. As a result of the alignment of the NLO maxima, in Fig. 6 we observe that the predictions based on the two scales that depend on the jet transverse energy, i.e. µ T,tot and µ T,jets , overlap almost perfectly, both at LO and NLO. A similarly good alignment is observed also between the other two scales, µ ttbbj and µ gbb , which do not depend on H T . Vice versa, the scales that do and do not depend on H T feature a non-negligible difference. In particular, the values of σ NLO at the maxima differ by about 10%. Such differences are most likely due to the fact that the dependence on H T , which is sensitive to NLO radiation, leads to a significant difference between the average scales in Born-like and real-emission contributions at NLO. Nevertheless, we observe that for all curves the position of the maximum coincides quite precisely with the intersection of the NLO and LO curves, which corresponds to K = 1. Moreover, the four K-factors coincide almost exactly in the whole ξ R range.
In summary, applying a rescaling that aligns dynamic scales based on the positions of the NLO maxima makes it possible to remove trivial differences related to the scale  normalisation and to highlight genuine differences related to their kinematic dependence. Since such alignment is in part already realised in the original scale choices (3.3),(3.5)-(3.7), in the following we will refrain from applying the small extra rescaling (4.13).

Fiducial cross sections
In this section we present detailed numerical results for integrated cross sections and scale uncertainties.
To highlight the quantitative importance of light-jet radiation emitted by the ttbb system, in Table 3 we present ttbb+jets cross sections with variable b-jet and light-jet multiplicities. Comparing the cross sections in the ttbbj and ttbb phase spaces, both available at NLO, we observe that the production rate for an extra light jet is around 50%, i.e. every second ttbb event involves a hard light jet with p T > 50 GeV. The ratio of the cross sections in the ttbbjj and ttbbj regions is around 40%, i.e. the emission of a second extra jet seems to be less abundant. However, one should keep in mind that this ratio is only LO accurate. The light-jet emission rates observed in the phase space with N min b = 1 are comparably large to the N min b = 2 case. For fixed N min j , cross sections with two b-jets are about a factor ten smaller wrt the corresponding cross sections with one b-jet. In general, LO scale uncertainties are very large, and grow by roughly 20% at each extra emission. Instead, scale uncertainties at NLO are drastically reduced, and in ttbbj production they are less pronounced than in ttbb production.
In the following we focus on LO and NLO predictions for ttbb+jet production in the ttbj and ttbbj phase-space regions. In Table 4 we compare cross sections and scale   Seven-point scale variations at LO are between around −45% and +90% for all scale choices, both in the ttbj and ttbbj regions. At NLO they are reduced around 20%, with significant differences depending on the scale choice and the number of b-jets. In the ttbbj (ttbj) phase space, the half-width of the scale-variation band is around 20% (25%) for the H T -independent scales and about 5% smaller for the H T -dependent ones.
In the last two columns of    Table 4 (15)%. We also note that the variation bands of the H T -independent scales cover the nominal predictions of the H T -dependent scales, but not vice versa. Based on this observations, we conclude that the somewhat larger seven-point variation of the H T -independent scales should be regarded as a more realistic estimate of scale uncertainties.
In Table 5 we present similar results based on the aligned scales (4.13), which correspond to Fig. 6. The main effect of the alignment is that the LO and NLO cross sections based on the two H T -independent scales become much closer to each other, while predictions based on the H T -dependent scales change in a less significant way. This is mainly due to the fact the original scales µ T,jets and µ T,tot are already very close to the corresponding aligned scales in (4.13). In any case, predictions based on the aligned scales are independent of the initial normalisation of the various scales.
After alignment, we still see significant differences between the predictions with H Tdependent and H T -independent scales. More precisely, due to the fact that the alignment is based on the NLO maximum of the ttbbj cross sections, the spread between K-factors in the ttbbj phase space goes down from 0.11 to 0.03. Vice versa, the K-factor difference in the ttbj phase space increases from 0.11 to 0.16. The alignment leads also to a slight reduction of NLO scale uncertainties, and the nominal predictions based on H Tindependent scales remain above the NLO bands of H T -dependent scales. Such differences between aligned NLO predictions in different phase-space regions should be regarded as genuine effects of the kinematic dependence of dynamic scales. Thus they play a largely complementary role wrt factor-two scale variations.

Sudakov effects
In this section we address the question of the safeness of the chosen transverse-momentum cut of 50 GeV with respect to higher-order Sudakov logarithms. To investigate such Sudakov effects, which appear in the region where the p T of the light jet, p T,j , becomes small, we relax the cut on p T,j and, in Fig. 7, we study the perturbative behaviour of the dσ/dp T,j distribution. In the left plot, this is done by keeping the usual b-jet cuts at p cut T = 50 GeV, while in the right plot this threshold is lowered to p cut T = 25 GeV. As is well known, the dσ/dp T,j distribution is logarithmically divergent at LO, while summing such logarithms to all orders in α S would cancel the divergence and lead to dσ/dp T → 0 at small p T . In the fixed-order NLO calculation at hand, this behaviour manifests itself through an increasingly strong shape difference between the LO and NLO distributions at small p T . For p cut T = 50 GeV, we find that at p T 20 GeV the NLO curve develops a Sudakov peak, below which NLO corrections start overcompensating the logarithmic growth of the LO distribution. In correspondence with the Sudakov peak, the NLO cross section is already less than half of the LO one, and below 15 GeV it rapidly dσ LO /dp T,j dσ NLO /dp T,j 1.45 0.881 0.699 1.09 0.754 0.639 Table 6. Comparison of the LO and NLO distributions in the leading-jet p T for p cut T = 50 GeV and 25 GeV. The results correspond to the ttbbj phase space with the cut p T > p cut T restricted to b-jets. falls into the unphysical regime of negative cross sections. This pathologic behaviour of the fixed-order NLO prediction is also reflected by the rapid inflation of NLO scale uncertainties below 40 GeV, while our choice of setting the light-jet p T cuts at 50 GeV guarantees good stability both for the NLO predictions and their uncertainties.
As can be seen in the right plot of Fig. 7, reducing the b-jet threshold to 25 GeV tends to lower the position of the Sudakov peak by 5 GeV or so. In this case, NLO predictions feature a good perturbative convergence down to 30-35 GeV. The effect of NLO corrections on the jet-p T distribution for selected values of p T is reported in Table 6.

Differential observables
In this section we study differential observables for pp → ttbbj at 13 TeV restricting ourselves to the ttbbj phase space. The main focus of our analysis is on the shapes of distributions and related uncertainties.

Distributions and shape uncertainties in the ttbbj phase space
In Figures 8-16 we analyse a series of differential distributions showing, for each observable, absolute and normalised distributions as well as six different ratio plots, which quantify the relative effects of seven-point variations and differences between the various dynamic scales. We restrict ourselves to the three dynamic scales (3.3),(3.5)-(3.6), since including or not the scale (3.7) does not change the overall picture of shape uncertainties. The format of the plots is described in the following and in the caption of (L2) A first ratio plot corresponding to the inverse K-factor, where scale variations are applied only in the numerator.
(L3) A second ratio plot that features the LO and NLO ratios, This ratio encodes differences between the dynamic scale µ R = µ gbb , defined in (3.5), and the default scale. Seven-point scale variations are applied in a correlated way to the numerator and the denominator. In this way, the main effect of factor-two variations, which amounts to a nearly constant normalisation shift, cancels out. As a result, the ratio (5.2) is mostly sensitive to effects that arise from the different kinematic dependence of the considered scales, and cannot be accounted for by factortwo variations of a single scale.
The right plot of each figure shows the following normalised distributions and ratios thereof.
(R1) The upper frame displays the LO and NLO normalised distributions, for the default scale µ R = µ ttbbj . The denominator corresponds to the integrated cross section in the ttbbj phase space, and seven-point variations in the numerator and denominator are correlated. In this way, distributions are always normalised to one, i.e. normalisation effects cancel out, and only shape corrections and uncertainties remain visible.
(R2) The first ratio plot shows the ratio of normalised distributions, (R3) The second ratio plot shows the ratios of normalised distributions at LO, for the three dynamic scales µ R = µ ttbbj , µ gbb , µ T,tot . This ratio highlights shape differences between those scales (with seven-point variations) and the nominal default scale.
(R4) The third ratio plot shows the same ratios as defined in (5.5), but at NLO,  3) based on the default scale, with correlated sevenpoint variations in the numerator and denominator. The first ratio plot (R2) displays the ratiô R (N)LO (µ ttbbj ), which is defined in (5.4) and highlights the relative shape distortions induced by NLO corrections and scale variations. The last two ratio plots on the right feature the ratios R (N)LO (µ R ) for µ R = µ ttbbj , µ gbb and µ T,tot at LO (R3) and NLO (R4). As defined in (5.5)-(5.6), such ratios quantify shape uncertainties associated with the kinematic dependence of the different dynamic scales. Fig. 8 presents the distribution in the p T of the leading light jet up to 400 GeV. The corrections to the shape of this distribution indicate excellent perturbative stability in the hard region above 150 GeV: The default scale yields a nearly constant K-factor around 1.65, and the scale-variation band is also quite stable at the ±20% level. In the region below 150 GeV, as already observed in Fig. 7, NLO effects start affecting the p T -shape with a correction of about 25% between 150 and 50 GeV. Such effects can be attributed to Sudakov logarithms, and estimating the missing higher-order corrections via naive exponentiation,  we expect residual shape uncertainties below 5% at NLO.
Comparing predictions based on the default scale and the other dynamic scales, in L3-L4 we observe normalisation differences at the level of 10-15%, which are compatible with the NLO scale-variation band in L2. These differences are very stable with respect to correlated factor-two scale variations as defined in (5.5): at LO such variations cancel almost exactly, and also the NLO bands in L3-L4 are suppressed at the level of 5% or less. Comparing normalised distributions with different dynamic scales in R3-R4, we see that LO shapes (and their seven-point variations) are almost identical, with only fewpercent differences between µ T,tot and the H T -independent scales. The nominal NLO predictions based on the various scales feature a similarly high level of agreement (see R4). However, similarly as in R2, factor-two variations lead to shape distortions at the 20% level. Such distortions shift the shape of the distributions in the region below 150 GeV, and are compensated by an opposite, but p T -independent shift in the hard region. In general, the suppression of shape effects at LO demonstrates the importance of NLO predictions for a more realistic assessment of shape uncertainties.
The non-negligible NLO shape effects observed in Fig. 8 are a specific feature of the jet-p T distribution in the vicinity of the cut, while other distributions that involve the leading light jet are typically more stable. This is illustrated in Figures 9-10 where we present the distributions in the pseudorapidity of the leading jet and in its ∆R separation with respect to the leading b-jet. For these observables, NLO corrections and uncertainties correspond to the ones of the  integrated cross section and depend only very weakly on the jet kinematics. In fact, as can be seen from the ratio plots R2-R4, the shape of such distributions turns out to be stable at the percent level with respect of seven-point variations and differences between dynamic scales.
In general, as found in Figures 8-10 and in various other observables not shown here, distributions in the leading light jet can be controlled with typical normalisation uncertainties of order 20% and shape uncertainties of order 10% or below.
In Figures 11-16 we present distributions in the top-quark and b-jet kinematics. For the transverse momentum of the harder top quark, shown in Fig. 11, we find that NLO corrections and scale variations are very stable, the only exception being a NLO shape correction of about 15% in the region below 50 GeV, where the cross section is strongly suppressed. For the p T of the softer top quark, shown in Fig. 12, NLO corrections feature a moderate, but more significant kinematic dependence. In particular, the K-factor goes down from about 1.5 in the bulk of the distribution to 1.2 in the tail, while seven-point scale variations lead to a similarly large shape distortion in the tail (see R2, R4). This behaviour is qualitatively quite similar to the Sudakov effects observed in the soft region of the jet-p T distribution in Fig. 8. It can be attributed to the fact that requiring two very hard top quarks restricts the available phase for additional radiation, confining the light jet into the soft region close to the 50 GeV threshold.
The distributions in the p T of the harder and softer b-jets, shown in Figures 13-14, feature a qualitatively very similar behaviour as the corresponding top-quark distributions.      In the case of the harder b-jet p T , NLO corrections and scale uncertainties depend rather weakly on p T (although more significantly than for the harder top quark), while the distribution in the p T of the softer b-jet features strong NLO effects, which are most likely due to Sudakov logarithms.
Finally, in Figures 15-16 we show the ∆R separation and the invariant-mass distribution of the b-jet pair. For these observables, as far as the default scale and the scale µ R = µ T,tot are concerned, NLO corrections and variations feature very little kinematic dependence, with percent-level shape differences. On the contrary, the dynamic scale µ R = µ gbb leads to a very different shape in the tail of the ∆R bb distribution, with deviations that reach −45% at LO and remain as important as −30% at NLO. A similar, although less dramatic trend is observed also in the tail of the invariant-mass distribution, which is clearly correlated to the tail of the ∆R bb distribution. These effects are most pronounced at ∆R bb > π, where the two b-jets are emitted in opposite hemispheres. In this region, the main mechanism of ttbb production via final-state g → bb splittings (see Fig. 3a) is strongly suppressed, and the leading role is played by topologies with initial-state g → bb splittings (see Fig. 3b in this paper and Fig.6 in [12]). The latter are maximally enhanced at E T,b m bb , and their characteristic virtualities of order E T,b are correctly reflected in the definition of the scales µ ttbbj and µ T,tot . Instead, the term m bb in (3.5) renders µ gbb unnaturally hard, leading to an unphysical suppression of the tails. It is clear that this behaviour cannot be regarded as a theoretical uncertainty, but should simply be taken as an indication that the scale µ gbb , which was designed to account for final-state g → bb splittings, is not applicable to initial-state g → bb splittings. On the contrary, the scales µ ttbbj and µ T,tot turn out to be well behaved for both kinds of splittings.

Recoil observables
As pointed out in the introduction, the accuracy of NLO Monte Carlo simulations of ttbb production plays a key role in ttH analyses. In this context, it was recently observed that the modelling of recoil effect by the parton shower may be a dominant source of uncertainty (see e.g. [50,51]). This is not surprising, given that every second ttbb event is accompanied by QCD radiation with p T > 50 GeV (see Table 3). In fact, away from the collinear regions, the recoil prescriptions used by parton showers can easily lead to unphysical momentum shifts of the order of 10 GeV and beyond. In the case of b-jets the effects of recoil mismodelling can be quite significant. In particular, shifts in the transverse momentum of the second b-jet can easily result in sizeable migration effects from the strongly populated region with N min b = 1 to the less populated N min b = 2 region. 10 In this context, the accurate description of QCD radiation provided by the calculation of 10 We have verified that in the ttbj region the second b-jet is typically slightly below the pT acceptance cut and is almost ten times softer with respect to the leading light jet. Thus, a small fraction of the QCD recoil is sufficient in order to shift the softer b-jet above the acceptance cut. More precisely, in the ttbj (ttbbj) phase space with standard cuts at 50 GeV the average transverse momenta of light jets and b-jets are pT,j 1 = 131 (137) GeV, p T,b 1 = 134 (166) GeV and p T,b 2 = 35 (86) GeV, while their average ratios are pT,j 1 /p T,b 1 = 1.34 (1.09) and pT,j 1 /p T,b 2 = 9.15 (1.83) . In the ttbj phase space, the quoted p T,b 2 and pT,j 1 /p T,b 2 averages have been evaluated including only events that involve a second resolved b-jet with p T,b 2 > 0 and |η b 2 | < 2.5. Figure 17. Sketch of the azimuthal angular correlation ∆φ rec,X between individual objects of the ttbb system and its recoil. See (5.7)-(5.8).
pp → ttbbj at NLO can be exploited as a benchmark to test the modelling of recoil effects in ttbb Monte Carlo simulations. With this motivation in mind, we study the azimuthal angular correlations [51] ∆φ rec,X = ∆φ ( p T,rec , p T,X ) (5.7) between the transverse momentum of the recoil, 8) and the various objects X of the ttbb system, i.e. the harder and softer top quarks (t 1 , t 2 ) and the harder and softer b-jets (b 1 , b 2 ), as well as the top-quark and the b-jet pairs. These angular observables, sketched in Fig. 17, reveal whether the respective object X absorbs a significant fraction of the QCD recoil through the presence (or absence) of peaks at ∆φ rec,X = ±π.
In Fig. 18 we present LO and NLO predictions for the azimuthal correlations between the recoil and the various top-quark and b-jet objects. For these observables we focus on the default scale, µ R = µ ttbbj , with seven-point variations. The absolute distributions in the upper frames indicate a very clear pattern: the recoil is preferentially absorbed by the harder top quark, and consequently also by the tt system, while the softer top quark and the b-jets feature only weak angular correlations with respect to the recoil. More precisely, in the case of the harder top, at ∆φ = ±π the cross section is almost five times larger as compared to the central region, while in the case of the harder (softer) b-jet this enhancement goes down to about 50% (20%). Thus it should be clear that naive shower models that distribute the recoil in a democratic way may lead to a significant mismodelling of the b-jet kinematics. Concerning the accuracy of NLO predictions in Fig. 18, we observe that all distributions are quite stable wrt to NLO corrections and scale variations. The most significant shape effects show up in the case of top-quark observables, where scale uncertainties can shift the level of the recoil peak by 15-20%, while for b-jets the flatness of the azimuthal correlations is remarkably stable with respect to higher-order effects.
These results demonstrate that fixed-order NLO predictions for pp → ttbbj can be used as a precision benchmark to validate the modelling of recoil effects in Monte Carlo simulations of ttbb production.  Figure 18. Distributions in the azimuthal angular separation ∆φ rec,X between individual objects X of the ttbb system and its recoil. See (5.7)-(5.8). The left column shows the angular correlations between the recoil and the top-quark objects X = t 1 , t 2 , t 1 t 2 , where t 1 t 2 denotes the top-pair system. Corresponding observables for b-jet objects, X = b 1 , b 2 , b 1 b 2 , are shown in the right column. Same setup and plots as in Fig. 7.

Tuning of QCD scale choice in ttbb production
In the literature on pp → ttbb at NLO, the usage of dynamic scales of type µ R = µ ttbb (3.4) has been advocated on the basis of the moderate size of the resulting NLO correction factor, K = σ NLO /σ LO . However, as pointed out in [12], the smallness of the observed Kfactor was largely due to the usage of a rather high LO value of α S as input for σ LO , while using the same α S in σ LO and σ NLO results in a correction factor as large as K 1.9 [12]. The lack of perturbative convergence, reflected by this large K-factor, may simply be the consequence of the fact that µ R = µ ttbb is a suboptimal choice. At the same time, it may also be the origin of the discrepancies between NLOPS simulations of ttbb production [52]. In fact, when matrix elements at NLO are matched to parton showers, the spectrum of the hardest QCD emission receives uncontrolled corrections of order (K − 1) = O(α S ). Such effects are formally beyond NLO, but for K 1 they can lead to sizeable distortions of the radiation spectrum [52].
In the light of these observations, and given the strong scale dependence of the ttbb K-factor, it is clear that a relatively mild reduction of the nominal scale would automatically lead to a smaller K-factor and, possibly, also to an improved behaviour of NLO matched simulations. However, the large ttbb K-factor may also be due to large higherorder effects that are not related to the choice of µ R . In this case, a reduction of the K-factor via µ R rescaling would only give a misleading impression of perturbative convergence without curing any problem. These considerations raise the question whether a reduction of the ttbb K-factor through a smaller choice of µ R may be supported through solid theoretical arguments. Generic considerations based on naturalness and perturbative convergence point towards a reduction of the standard ttbb scale choice by a factor 1/2 to 1/3 [52]. However, only the knowledge of the next perturbative order can shed full light on the goodness of a scale choice, i.e. on its effectiveness in capturing the dominant higher-order effects. In the case at hand, the ttbb scale choice could be tuned based on the requirement that i.e. by optimising the choice of the scales µ opt R,F in such a way that NLO ttbb predictions match NNLO ones. 11 However, the required NNLO calculation is completely out of reach. Nonetheless, the NLO corrections to pp → ttbbj presented in this paper represent one of the building blocks of ttbb production at NNLO, and as such they can provide useful insights on how to improve the ttbb scale choice. The idea is that the condition (6.1) can be imposed at the level of the jet-radiation spectrum by requiring With other words, the scale choice can be tuned in such a way that the tree-level description of the jet-p T spectrum that results from the ttbb NLO calculation matches the more precise prediction of the ttbbj NLO calculation. Contrary to (6.1), this procedure cannot guarantee the correct description of higher-order effects at the level of the inclusive ttbb cross section. 12 Nonetheless it is attractive for at least two reasons. First, tuned ttbb NLO predictions will guarantee a much more accurate description of the jet-p T spectrum, which is known to play a critical role in Monte Carlo simulations. Second, the shape of the jet-p T spectrum can be used to judge the quality of the matching procedure (6.2), and the general consistency of the procedure can be validated by comparing various other jet observables.
The results of this tuning procedure are presented in Fig. 19, where we show the distribution in the p T of the hardest light jet, and in the invariant masses of the systems formed by the hardest light jet in combination with the leading or the subleading b-jet.
The tuning is carried out through a constant rescaling of the standard ttbb scale choice, such as to match NLO predictions for the integrated ttbbj cross section based on the default scale µ ttbbj . To be conservative, we have compared two possible ways of tuning the ttbb scale. In the first approach, the rescaled ttbb NLO predictions are matched to nominal ttbbj NLO predictions, whereas in the second approach the tuning is done by matching the average values of the respective seven-point variation bands. The outcome of these two matching prescriptions is shown in the left and right columns of Fig. 19. Matching nominal predictions leads to a reduction of the default ttbb scale by a factor 13 κ = 1/1.6, whereas matching the scale-variation bands in a symmetric way requires a significantly smaller rescaling, κ = 1/1.14. This large difference is mainly due to the strong asymmetry of the factor-two variation band of the tree-level prediction, i.e. pp → ttbb at NLO. In this respect, we note that such asymmetry is mainly due to the logarithmic nature of the scale dependence (4.2). Thus the asymmetry of the LO band would largely disappear on logarithmic scale, and the prescriptions based on the central scale and the average of the bands would be significantly closer to each other. For all considered jet observables we find that both tuning scenarios lead to a very good agreement, not only in the normalisation, but also at the level of shapes. The findings of this analysis support a reduction of the standard ttbb scale (6.3) by up to a factor κ ∼ 1/1.6. In the ttbb (ttb) phase space, κ = 1/1.6 corresponds to a reduction of the ttbb K-factor from 1.80 (1.92) to 1.51 (1.62) and an increase of the nominal ttbb cross section by 18% (21%). 12 We note that this approach does not improve the precision of the integrated ttbb cross sections. Its goal is only to optimise the choice of the central scale. 13 We have checked that keeping µF = HT/2 fixed and tuning only µR would require a rescaling factor κ = 1/1.76.  Figure 19. Distributions in the p T of the leading jet and the mass of light-jet-b-jet systems in the ttbbj phase space. Comparison of NLO ttbbj predictions with default scale choice, (µ R , µ F ) = (µ ttbbj , H T /2), to NLO ttbb predictions with (µ R , µ F ) = (κ µ ttbb , κ H T /2). In the left plots, the reference curves for the matching procedure (solid, labelled NLO) correspond to the above central scales, and the applied rescaling factor is κ = 1/1.6. In the right plots, the reference curves (solid, labelled NLO) are the average values of the scale-variation bands, and κ = 1/1.14. The blue dashed curves indicate the position of the NLO ttbb reference prediction before tuning, while all other NLO ttbb predictions and scale-variation bands correspond to the tuned scales.

Summary
Measurements of ttH(bb) production at the LHC require very accurate theoretical simulations of the irreducible ttbb background. To address the dominant sources of systematic uncertainties, which stem from the modelling of QCD radiation in ttbb events, we have presented a calculation of ttbb production in association with one extra jet at NLO QCD.
To carry out this non-trivial calculation we used OpenLoops 2 in combination with the Sherpa and Munich Monte Carlo frameworks. Technically, the calculation of the required 2 → 5 one-loop amplitudes has confirmed that the new algorithms implemented in OpenLoops 2 can tackle multi-particle and multi-scale problems with very high CPU efficiency and numerical stability.
We have discussed pp → ttbbj at the 13 TeV LHC with emphasis on the effects of NLO corrections and scale uncertainties. To this end, we have studied conventional factortwo rescalings, as well as variations of the kinematic dependence of dynamic scales. In order to disentangle normalisation and shape uncertainties in a transparent way, we have proposed to compare dynamic scales upon alignment of the NLO maxima of the respective scale-variation curves.
In general, the typical level of scale uncertainties in pp → ttbbj at NLO is 20-30% for integrated cross sections and below 10% in the shapes of distributions. The calculation at hand can thus be used as a precision benchmark to validate the modelling of QCD radiation in Monte Carlo generators of ttbb production. With this motivation in mind, we have presented NLO predictions for various azimuthal correlations that provide a transparent picture of the effects of the recoil of QCD radiation on the different objects of the ttbb system.
Finally, we have discussed the issue of the large NLO K-factor observed in inclusive NLO calculations of ttbb production, and we have addressed the question of whether it is justified to reduce this K-factor through ad-hoc scale choices. In this respect we have argued that the NLO corrections to pp → ttbbj entail information on pp → ttbb beyond NLO, which can be exploited to identify an optimised scale choice. Specifically, we have proposed the idea of adjusting the nominal ttbb scale choice such as to match the jet emission rate predicted by pp → ttbbj at NLO. This improved scale choice leads to a reduction of the ttbb K-factor, and is also expected to attenuate theoretical uncertainties in the context of NLO matching to parton showers.