Production and decay of the Higgs boson in association with top quarks

We report on the calculation of the next-lo-leading order QCD corrections to Higgs boson production and decay in association with top quarks. We consider leptonic decays of top quarks leading to the hadronic process pp → e+νeμ−ν¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\nu} $$\end{document}μbb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}H(H → X) at the LHC with s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV. All resonant as well as non-resonant Feynman diagrams, interferences and off-shell effects are included for the top quark and W gauge boson. Decays of the Higgs boson, on the other hand, are included in the narrow-width approximation. Specifically, we consider Higgs boson decays into bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}, τ+τ−, γγ and e+e−e+e−. Numerical results are given at the integrated and differential fiducial level for various factorisation and renormalisation scale choices and different PDF sets. We study the main theoretical uncertainties that are associated with neglected higher order terms in the perturbative expansion and with different parametrisations of the PDFs. Furthermore, we examine the size of the off-shell effects by an explicit comparison to the calculation in the full narrow-width approximation. Finally, the impact of the contributions induced by the bottom-quark parton density is investigated.


Introduction
After the discovery of the Higgs boson in 2012 by the ATLAS and CMS collaborations [1,2], one of the main goals of the Large Hadron Collider (LHC) is the precise determination of its properties and its couplings to other Standard Model (SM) particles. The Higgs boson couples to fermions via the Yukawa interaction, which is proportional to the mass of fermions and therefore of particular interest for the top quark, the heaviest fermion in the SM. Within the SM, the top quark Yukawa coupling (Y t ) is close to unity. Its measurement can not only serve as an important test of the SM but also can be used to constrain various models beyond the SM (BSM). Indeed, in many new physics scenarios Y t is expected to differ from the SM value. Higgs production at the LHC is dominated by gluon fusion, gg → H. Since this process is loop induced, it only provides an indirect way to access Y t , and it can be affected by heavy BSM particles. On the other hand, Higgs production in association with top quarks, pp → ttH, gives a window to directly probe this coupling already at tree level. This production mechanism contributes only about 1% [3] to the total Higgs boson production cross section and was observed in 2018 by the ATLAS and CMS collaborations [4,5].
On the other hand, the Higgs boson can only be observed in the detector by its decay products. The first observation of ttH production in a single decay channel was reported recently by the ATLAS and CMS collaborations [6,7] through Higgs decay into two photons JHEP02(2022)196 H → γγ. Even though this decay channel is loop induced and has a very small branching ratio of about 0.2% [3], it is experimentally very well accessible due to the small background and good mass resolution for M γγ . Measurements of the Higgs boson production rate in association with top quarks in final states with electrons, muons and tau leptons have been carried out very recently by CMS [8]. This analysis aimed at events that contain H → W + W − , H → τ + τ − or H → ZZ decays with the corresponding branching ratios of about 21%, 6% and 3% respectively. Further searches are reported by ATLAS and CMS for ttH with a Higgs boson decay into a bottom quark pair H → bb [9,10], which has the largest branching ratio of about 58% [3]. From the experimental point of view, however, this decay channel is the most challenging one, because of the complicated Higgs boson reconstruction from the so-called combinatorial background caused by a multitude of bjets coming both from the two top quarks and the Higgs boson. Furthermore, this decay channel suffers from large reducible ttjj background and irreducible QCD and Z-peak backgrounds, ttbb and ttZ(Z → bb) respectively.
In this article, we report on the computation of the NLO QCD corrections to offshell ttH production in the di-lepton top quark decay channel with Higgs boson decays in the NWA. In a first step, we compute independently the NLO QCD corrections to pp → e + ν e µ −ν µ bb H + X at order O(α 3 s α 5 ) at the LHC with √ s = 13 TeV and perform a comparison with the results presented in ref. [33] at the integrated and differential fiducial level. In a next step, we extend the analysis of ttH production and investigate the dependence of our results upon variation of renormalisation and factorisation scales and parton distribution functions in the quest for an accurate estimate of the theoretical uncertainties. Additionally, we explore a few possibilities for a dynamical scale choice with the goal of stabilizing the perturbative convergence of the differential cross sections in the high p T tails. Afterwards, we study the off-shell effects of the top quarks at the integrated and differential level by a comparison with the calculations in the NWA and NWA with leading order (LO) top quark decays (NWA LOdecay ). Having NWA LOdecay results at our disposal, we can additionally estimate the size of the NLO QCD corrections to top-quark decays. We close the discussion on ttH production by investigating the b-initiated contributions that are often neglected in similar calculations with off-shell top quarks because of their assumed numerical insignificance. At last, we include decays of the Higgs boson in the NWA while we keep all off-shell effects of the top quark and W gauge boson. Even though the non-factorisable corrections vanish in the limit Γ H /m H → 0, that characterises the JHEP02(2022)196 NWA, such a mixed approach is very well justified when we observe that the following hierarchy is satisfied 2.6% > 0.8% 0.003% .
In our paper, various decay channels of the Higgs boson are added. Thus, we can analyse fully realistic final states. We study the scale dependence in the decays of the Higgs boson and compare it to the one in the Higgs boson production. Additionally, various differential distributions of Higgs boson decay products are investigated in more detail. The paper is organised as follows. In section 2, we briefly describe the Helac-Nlo framework used for the ttH process, and the computation of Higgs boson decays, as well as the various cross-checks performed in this calculation. The theoretical setup, consisting of input parameters, event selection, renormalisation and factorisation scale choices, is described in section 3. In section 4 we perform a comparison with the results of ref. [33]. Phenomenological results for ttH production at the integrated and differential fiducial level are discussed in detail in section 5 and section 6. In section 7 we examine the off-shell effects of the top quarks. The contribution from bottom quarks in the initial states and the modifications required in the jet algorithm are discussed in section 8. The production and decay of the Higgs boson are combined in section 9. We consider Higgs boson decays into bb, τ + τ − , γγ and e + e − e + e − , and present results at the integrated and differential fiducial level. Finally, our results are summarised in section 10.

Description of the calculation
In a first step, we study Higgs boson production in association with top quarks while keeping the Higgs boson stable (Γ H = 0). The complex-mass scheme [36][37][38][39] is used to include all off-shell effects of the top quarks and massive gauge bosons (W ± , Z) in a gauge invariant way. In particular, we consider the following process pp → e + ν e µ −ν µ bb H + X the Dyson-Schwinger off-shell iterative algorithm [40][41][42] instead of the traditional Feynman diagrams approach. This iterative algorithm is implemented in Helac-Dipoles [43] based on Helac-Phegas [44][45][46]. The phase-space integration is performed and optimised with the help of Parni [47] and Kaleu [48]. The virtual corrections are computed with Helac-1Loop [49], where the one-loop matrix elements are reduced to scalar integrals at the integrand level using the OPP reduction technique [50] as implemented in Cut-Tools [51]. The scalar integrals are obtained from the library OneLOop [52]. Since the computation of loop diagrams is very time consuming, the number of evaluated loop diagrams is reduced by using only unweighted Born events for the calculation of the virtual corrections as described in ref. [53]. The finite part and -poles are cross-checked with the module MadLoop [54] in MadGraph [55] for a few phase-space points. In addition, the cancellation of the -poles between the virtual corrections and the I-operator implemented in Helac-Dipoles is checked. For the real corrections we have to take into account the same sub-processes as above but with an additional gluon in the final state. We have 1162 Feynman diagrams for the gluon induced gg → e + ν e µ −ν µ bb Hg and 476 for the quark induced qq → e + ν e µ −ν µ bb Hg sub-process. At NLO, we encounter additional sub-processes, namely gq (gq) → e + ν e µ −ν µ bb Hq (q), which are connected by crossing with the quark induced ones. The calculation of the real corrections is performed via subtraction methods implemented in Helac-Dipoles and is identical from the QCD point of view with the calculations in refs. [56][57][58]. In particular, two independent schemes are implemented, the Catani-Seymour dipole formalism [59,60] and the Nagy-Soper subtraction scheme [61]. In addition, a phase-space restriction on the subtraction terms is implemented for both schemes as described in refs. [62][63][64][65]. This restriction is used to reduce the number of required subtraction terms and as a cross-check, since the final real emission results should not depend on the restriction. For the full off-shell calculation the Nagy-Soper subtraction JHEP02(2022)196 scheme is used, while for the NWA case the Catani-Seymour dipole formalism is employed instead. The calculation in the NWA is also performed in the Helac-Nlo framework and details about the implementation are given in ref. [66]. All parts of the off-shell calculation are stored in modified Les Houches Event files (LHEFs) [67] as partially unweighted events [68] with additional information based on ref. [69]. This allows us to create differential distributions and obtain results for different PDF sets or other factorisation and renormalisation scale settings without rerunning the entire calculation from scratch.
In the last part, we combine the production and the decay of the Higgs boson with a tt pair. We reuse the LHEFs to model the Higgs boson decays in the NWA. Consequently, we are able to include the production and decays of the Higgs boson in association with top quarks using the NWA for the Higgs boson at NLO in QCD, while keeping full offshell effects for top quarks and for all other unstable particles. This approach is very well justified as we argued before. As we consider decays of a scalar particle, no spin information is transferred from the production to the decay stage which further warrants the use of the LHEFs. The phase-space integration for Higgs boson decay products is performed using our own implementation of the algorithm outlined in ref. [70], called Rambo (on diet). In particular, we generate unweighted events in the rest frame of the Higgs boson to model its decays. The implementation is cross-checked at LO with the NWA implementation in Helac-Dipoles. All parts of the calculation with Born-like kinematics work exactly in the same way. In this paper we consider the following decay channels of the Higgs boson In the last case, the Higgs boson decays into two off-shell Z-bosons, which further decay into two e + e − pairs [30]. Since we are interested in the e + e − e + e − events in the fiducial phase-space regions near the Higgs boson mass peak, this approach is very well justified. For the 4 events well above the Higgs boson mass peak, i.e. when M 4 m H , sizable Higgs boson continuum interference effects might occur that need to be considered as well, see e.g. refs. [31,32]. We do not include leptonic τ ± and W ± -gauge boson decays, even though the two decay chains seem to be straightforward in the NWA approach. The reason is rather technical, as we simply tend to avoid the final states were additional neutrinos can appear. This makes it easier to directly compare with the results presented in ref. [33] where a missing transverse momentum cut has been applied to the two neutrinos from the top quark decays already at the Higgs boson production stage. Of course, had we generated the LHEFs without this cut in the first place, both H → τ + τ − → + ν −ν ν τντ and H → W + W − → + ν −ν with ± = e ± , µ ± could have been added as well. To keep all obtained events and to avoid generation of new LHEFs, we decided to focus only on the final states without neutrinos. It should be noted that generating such a complicated 2 → 7 final state is very time consuming and requires a dedicated computer cluster. The JHEP02(2022)196 differential cross section for the pp → ttH → ttX process, where tt symbolises the following 2 → 6 process pp → e + ν e µ −ν µ bb and X stands for X = bb , τ + τ − , γγ , e + e − e + e − , can be written as In the second line the Higgs boson production cross section and decay rate are expanded up to NLO in QCD, while Γ H is kept fixed. Specifically, dσ 0 ttH , dΓ 0 H→X represent LO contributions as calculated with the NLO input parameters and dσ 1 ttH , dΓ 1 H→X describe the NLO corrections. The last term is only present in the case of the Higgs boson decay into a bb pair. 1 For H → bb also the I-operator and the virtual corrections are calculated using LHEFs. We use analytical results for the virtual corrections and the I-operator where the latter one is explicitly cross-checked with the I-operator implemented in Helac-Dipoles. The Catani-Seymour dipole formalism with phase-space restriction on the dipoles, as implemented in Helac-Dipoles, is used for the calculation of the real corrections in this decay chain. For the H → γγ decay, NLO QCD corrections are neglected as these would require two-loop integrals. Furthermore, it is well know that the NLO corrections relative to the lowest order are small in general as shown in ref. [71]. As we previously mentioned, the five flavour scheme is used. Thus, the bottom quark mass is still neglected throughout the calculation. This allows us to resum large initial state collinear logarithms into PDFs. In order to keep a non-zero coupling between the Higgs boson and bottom quarks, we keep Y b non-zero in the decay stage and perform the renormalisation in the MS scheme. Moreover, we distinguish between the renormalisation scale in the production µ R and decay µ R, dec . The decay width into a bottom quark pair at NLO is checked with the result in ref. [72] by setting µ R, dec = m H . In H → τ + τ − decays, we do not neglect the mass of the τ lepton in the matrix element or the phase space. For the loop-induced decay into two photons we use the matrix element at LO as given in refs. [73,74]. This matrix element comprises two types of contributions: from W gauge bosons and from massive fermions. In the latter case we consider massive fermion loops from top quarks, bottom quarks and τ leptons. However, the contribution from the latter two is negligible.

LHC setup
We consider the pp → e + ν e µ −ν µ bb H + X process at the LHC with √ s = 13 TeV. We use the same input parameters as in ref. [33], as briefly summarised in the following. The 1 We note that terms like dσ 1 ttH dΓ 1 H→X are not included in this defintion since such terms would be of the order of O(α 4 s ) and are part of NNLO corrections. With this definition, the total NLO cross section is not equal to the NLO production cross section times the NLO branching ratio. On the other hand, this relation can be obtained by simply rescaling the terms in eq. (2.2) by appropriate factors. However, such a definition might lead to underestimated theoretical uncertainties and even affect the scale dependence of the production process as outlined for example in ref. [87] (see appendix B) for the pp → V H + X → + − bb process, where V = W ± , Z.

JHEP02(2022)196
number of active flavours is set to N f = 5 and the contribution from bottom quarks in the initial state is neglected at first. We will return to this contribution later in this article. The CT10NLO [75] PDF set is used at LO and NLO QCD for the comparison, while newer PDF sets are employed for the rest of the paper. Specifically, as our default PDF set we use NNPDF3.1 NLO (LO) [76] with α s (m Z ) = 0.118. To quantify the differences between various PDF sets, that originate among others from the choice of the data used and the theoretical assumptions made for the global fit, we also provide numerical results at NLO QCD for CT18 [77] and MSHT20 [78]. The running of the strong coupling constant α s , with two-loop (one-loop) accuracy at NLO (LO) is provided by the LHAPDF interface [79]. The electromagnetic coupling α, is derived from the Fermi constant in the G µ −scheme The values for the masses and widths are given by where m OS V and Γ OS V with V = W, Z stand for the measured on-shell values (OS) for the masses and widths of the W and Z boson. They are converted into pole values according to ref. [80] The latter are used as input parameters for this calculation. The masses of all other quarks and leptons are set to zero. A Higgs boson mass of m H = 126 GeV is used to facilitate the comparison with results from ref. [33]. Even though the current value is closer to m H = 125 GeV, our conclusions do not depend on the specific value used. 2 Furthermore, the top quark width is derived from the equations in ref. [39] using NLO corrections as calculated in ref. [81]. The top quark width is treated as a fixed parameter throughout this work. Its value corresponds to a fixed scale µ = m t , that characterises the top quark decays. The α s (m t ) parameter is independent of α s (µ 0 ) that enters the matrix element calculations as well as PDFs, since the latter describes the dynamics of the whole process. The value of α s (m t ) is taken from the CT10NLO PDF set. In this way we obtain the following values for the top quark width in the full off-shell case whereas for the NWA case they are given by (3.5) 2 We have checked that for the integrated LO cross section changing the Higgs boson mass to mH = 125 GeV amounts to 2%. Indeed, for the NNPDF3.1 PDF set using µ0 = HT /2 one would obtain σ LO = 2.2618(4) fb instead of σ LO = 2.2130(2) fb.

JHEP02(2022)196
Final state partons with pseudo-rapidity |η| < 5 are clustered into jets using the IR-safe anti-k T jet algorithm [82] with the jet-resolution parameter R = 0.4. We require two charged leptons, two b-jets, missing transverse momentum and one stable Higgs boson. All final states have to pass the following experimental cuts where = e, µ. The factorisation, µ F , and renormalisation, µ R , scale are set to a common value µ F = µ R = µ 0 . For µ 0 we use the two scale settings defined in ref. [33] and In addition, we use a third scale setting where b 1 and b 2 are the two bottom-flavoured jets with highest transverse momentum. 3 This specification becomes relevant only when the contributions induced by the bottomquark parton density are included as well. In that case, we can have up to three b-jets in the final state. On the other hand, p T, miss is the total missing transverse momentum from escaping neutrinos defined by The additional scale setting is chosen such that information about the underlying resonant nature of the process is not used at all. It seems a more natural choice for the process at hand where full off-shell effects are included, thus, also Feynman diagrams with single-and non-resonant top-quark contributions are present at the matrix element level. As differential cross sections extend up to large energy scales, single-and non-resonant contributions might become significant, sometimes even dominant, in the specific phase-space regions, see e.g. ref. [66]. Thus, a H T based scale setting is better motivated, especially for dimenionful observables, than any scale choice which assumes the double-resonant nature of the process for each phase-space point. Furthermore, µ 0 = H T /2 does not require the reconstruction of the tt pair which is not free from ambiguities. Nevertheless, for the integrated fiducial cross section, which is dominated by the double resonant top-quark contributions, µ 0 = H T /2 leads to similar results as µ 0 = µ dyn and µ 0 = µ fix . Theoretical uncertainties arising from missing higher orders are estimated using a 7-point scale variation, in which the factorisation and renormalisation scales are varied independently in the range

JHEP02(2022)196
In practice, this amounts to considering the following pairs (1,2) , (1, 1), (1, 0.5), (2, 2), (0.5, 0.5) . (3.12) As it is usually done, we search for the minimum and maximum of the resulting cross sections. For the PDF error, we use the prescription of the NNPDF3.1 group to provide the 68% confidence level (C.L.) PDF uncertainties. Specifically, the NNPDF3.1 group uses a Monte Carlo sampling method in conjunction with neural networks, and the PDF uncertainties are obtained with the help of the replica method. On the other hand, a Hessian representation is used for CT18 and MSHT20, where the PDF uncertainties of CT18 have to be re-scaled to obtain the 68% C.L. since they are originally provided at 90% C.L. In the following, we summarise the additional parameters required for the decay of the Higgs boson. We use the masses given in ref. [3] for the τ lepton and the bottom quark in the on-shell and MS-scheme The Higgs width is calculated using the parameters described above and the program Hdecay [83,84], which leads to at NLO, where we have used m b (m b ) as initial value. We use the same kinematical cuts as already discussed for all Higgs boson decay products, including τ leptons. In addition, we require the following cuts for the two photons and apply Frixione's photon-isolation condition as defined in ref. [86] before the jet algorithm. In particular, for each parton i we evaluate R γi between this parton and the photon. We reject the event unless the following condition is fulfilled

JHEP02(2022)196
where we set γ = 1, n = 1 and R γ, j = 0.3. Furthermore, E T, i is the transverse energy of the parton i and E T,γ is the transverse energy of the photon. Jets reconstructed inside the cone size R γ, j are not subject to additional selection criteria. We note that the condition has to be satisfied for both photons from the H → γγ decay. For Higgs boson decays we still use the following scale setting µ 0 = H T /2 for µ R and µ F . However, now the transverse momentum of the Higgs boson needs to be reconstructed from its decay products. In the case of H → γγ, τ + τ − , e + e − e + e − , the Higgs boson can be reconstructed exactly since we have no corrections in the decays. Only for the H → bb decay we need a Higgs reconstruction prescription. In particular, we minimise the following where i, j run over all b-jets that are present, to find the two that originate from the Higgs boson decay. Morover, in eq. (3.17) M b i b j stands for the reconstructed Higgs mass and m H for the Higgs mass. The remaining two b-jets are assigned to the top quark decays. We choose to set and vary the renormalisation scales independently for the production and decay parts. We evaluate the cross section for a total of 21 different scale settings as described in ref. [87]. The latter are obtained from all possible combinations of where ξ ∈ [0.5, 1, 2] with the additional constraint following the conventional 7-point scale variation for the Higgs boson production process.

Numerical checks
To verify our results for the pp → e + ν e µ −ν µ bb H process, we perform a comparison of our predictions with the results from ref. [33] at the integrated and differential fiducial level. Since this calculation was performed with a completely independent computational framework, agreement of the results provides a strong indication of the correctness of the two calculations. In table 1, the integrated fiducial cross section and the K-factor, defined as K = σ NLO /σ LO , is shown for both calculations at LO and NLO QCD. Theoretical predictions are provided for µ 0 = µ dyn and µ 0 = µ fix for the CT10NLO PDF set. Also reported are theoretical uncertainties as obtained by 7-point scale variation together with Monte Carlo errors. We find good agreement between the two calculations at LO and NLO QCD for both scale settings. In particular, the K-factor and the scale uncertainties are also in agreement within the statistical uncertainties. The scale uncertainty comparison is an indirect comparison of different scale choices implied by the scale variation.
In figure 2, we perform a similar comparison at the differential fiducial level for the invariant mass of the bb system, M b 1 b 2 , the transverse momentum of the positron, p T,e + and the difference in azimuthal angle between the two charged leptons,    Results are evaluated using µ 0 = µ dyn and the CT10NLO PDF set. Also given is a comparison between Helac-Nlo and the results from ref. [33] (indicated by DF). Lower panels display the difference between both calculations in standard deviations (denoted as σ-dev. TeV. Results are evaluated using µ 0 = µ dyn and µ 0 = µ fix for the CT10NLO PDF set. Also given is a comparison between Helac-Nlo and the results from ref. [33] (indicated by DF).
at LO (left) and NLO (right). We employ µ 0 = µ dyn and the CT10NLO PDF set. For each observable, bin-by-bin, we display in the lower panels the difference between the two calculations in standard deviations (denoted as σ-dev.). The numerical values of ref. [33] were obtained using the online tool WebPlotdigitizer [88]. This tool can be used for the semi-automatic extraction of numerical values from images of data visualisations which makes an extraction of a large number of data points fairly easy. Potential uncertainties due to the extraction are neglected and the same Monte Carlo errors as in our calculation are assumed. Again we find overall good agreement between the two calculations at LO and NLO.

Integrated fiducial cross sections
In this section, we first analyze the two main sources of theoretical uncertainties resulting from the scale dependence and the PDF uncertainties at the integrated fiducial level. We extend the discussion of the scale dependence of ref. [33] by an additional dynamical scale setting µ 0 = H T /2 and investigate the differences between a fixed and dynamical scale choice. Second, we present results for the integrated cross section for more recent PDF sets and discuss the internal PDF uncertainties and the differences between the results for various PDF sets.
In table 2, we show the integrated cross section at LO and NLO QCD together with the corresponding K-factor for all three scale choices for the NNPDF3.1 PDF set. The deviations between the three scales are up to 5% at LO and less than 1% at NLO. The deviations are well within the scale uncertainties at the specific order in perturbation theory. Estimating the theoretical uncertainties by using the maximum of the lower and upper bounds from the scale variation, we find that they are very similar in size for all three scale settings. The scale dependence reduces drastically from about 30% at LO JHEP02(2022)196 to 5% at NLO, which corresponds to a reduction by a factor of 6. 4 Had we varied the factorisation and the renormalisation scale simultaneously instead, the maximum of the lower and upper bound of the scale variation would not change. The K-factor ranges from 1.18 to 1.23. The largest value is observed for µ 0 = H T /2 and is caused by the smaller value of the integrated cross section at LO. In summary, all three scales lead to very similar results at the integrated level and, in particular, we do not observe any differences between a fixed and a dynamical scale setting, which is not surprising as the latter should mainly play an important role at the differential level.
In figure 3, the integrated fiducial cross section at LO and NLO in QCD is shown as a function of the factorisation/renormalisation scale for all three scale choices for the NNPDF3.1 PDF set. In the first plot, the integrated cross section is shown for the simultaneous variation of the factorisation and renormalisation scale. The remaining three plots display the simultaneous variation and the variation of one scale, factorisation or renormalisation, while the other scale is kept fixed at the corresponding central value, µ 0 . When both scales are varied at the same time, we observe that for µ 0 = µ fix and µ 0 = µ dyn , the scale dependence is almost identical over the entire range both at LO and NLO. For the fixed scale, we observe that the integrated cross section becomes negative for small values at NLO. However, the choice of the PDF set plays a role since this behaviour is not observed for CT10NLO. For µ 0 = H T /2, our default scale choice, we find that at LO this scale setting leads to smaller values of the integrated cross section, where the deviations to the other two scales decrease for larger values. At NLO, all three scales are very similar in the range from 0.5 to 2 where the theoretical uncertainties arising from the scale dependence are actually evaluated. We can further notice that in this range, the dependence on the factorisation scale is constant. Thus, the scale variation is driven solely by the changes in µ R . Consequently, the 3-point and 7-point scale variation would lead to essentially the same results.
In figure 4, the scale dependence is shown for the normalised integrated cross section at LO and NLO in the µ R − µ F plane as a contour plot for µ 0 = H T /2 and µ 0 = µ fix . As the JHEP02(2022)196 [fb] scale dependence of µ 0 = µ fix and µ 0 = µ dyn is almost identical, the latter one is not shown again. This graphical visualisation allows a more detailed look at the scale dependence and, in particular, at the interplay between the changes in µ R and µ F . At LO and NLO, we again observe that changes in the renormalisation scale dominate the scale dependence. The latter is reduced significantly when moving from LO to NLO. In addition, at LO, we find only minor differences in the scale dependence between the two scales, while at NLO, larger effects can be found, especially for small values of the scales.
In addition to the scale uncertainties, it is important to study the second main source of theoretical uncertainties, that comes from the parameterization of PDF sets. We should note that there has been a lot of progress in determining the PDFs in the last few years. For example, various new NNLO processes have recently been included in PDF fits. In addition, all groups included large sets of LHC data, which resulted in significant changes to the global PDF fits. With increased statistical precision of many measurements, challenges JHEP02(2022)196  have arisen in fitting some of these new data sets. Moreover, differences in the fitting methodology, e.g. due to the treatment of correlated systematics, have a significant impact on the PDFs. Due to these advances, differences might emerge in both the central values and the uncertainties of the PDFs extracted by different groups. Consequently, it is vital to provide results for various PDF sets and study their internal uncertainties. In this way, it is possible to show the size of these differences, if any, and determine whether they are significant. In table 3 we provide the integrated fiducial cross section at NLO in QCD for NNPDF3.1, CT18 and MSHT20. The results are given for µ 0 = H T /2. Also shown are the corresponding scale and PDF uncertainties. We observe that the NNPDF3.1 PDF set yields the smallest PDF uncertainties of about 1.1% compared to 2.4% for CT18 and 1.9% for MSHT20. The PDF uncertainties are about a factor of 2-4 smaller than the corresponding scale uncertainties. Thus, scale uncertainty is still the dominant source of theoretical error. Furthermore, all three PDF sets lead to very similar predictions at the integrated level. Indeed, when comparing σ NLO   find the largest differences between NNPDF3.1 and CT18 of about 2.5%. We conclude that the PDF's internal uncertainties are similar in size to the differences between the three PDF sets.
In figure 5, the integrated cross section is shown at NLO QCD for µ 0 = H T /2 for the three previously used PDF sets. In addition, we provide theoretical results for their former versions, NNPDF3.0 [89], CT14 [90] and MMHT14 [91]. For each PDF set internal uncertainties are shown in black. For NNPDF3.1 the scale uncertainties are also shown in blue as reference. We find for all three groups that the former and current PDF sets agree well within their corresponding PDF uncertainties. The latter are reduced by about 1% for CT18 and by about 0.5% for MSHT20 and NNPDF3.1. We find overall good agreement between the different PDF sets, with the largest deviations between CT18 and NNPDF3.0

JHEP02(2022)196
of about 1.6σ or 4.4%, which is about the same size as the scale uncertainties and thus not negligible. The largest PDF uncertainties are found for CT14 of about 3.5%, which is still smaller but of comparable size to the scale uncertainties.

Differential fiducial cross sections and PDF uncertainties
We start this section with a comparison between the different scale settings at the differential level. We then continue the discussion of NLO QCD corrections and perturbative stability. At the integrated fiducial level, we have concluded that the scale choice has only a minor impact on the cross section. At the differential level, however, a fixed scale choice cannot handle the dynamics of the process in certain phase-space regions, leading to perturbative instabilities. It is well known that this problem can be solved with a judicious choice of a dynamical scale. This aspect is discussed below for various observables. We end this section with a comparison of the two main theoretical uncertainties.
We start in figure 6 with the differential distributions for the transverse momentum of the hardest b-jet (p T, b 1 ), the scalar sum of the transverse momenta of the two b-jets , the invariant mass of the first and second hardest b-jet (M b 1 b 2 ) and the invariant mass of the two charged leptons (M e + µ − ). NLO QCD results are shown for all three scale choices, µ 0 = µ fix , µ 0 = µ dyn and µ 0 = H T /2. Moreover, the NNPDF3.1 PDF set is employed. Scale uncertainties are displayed as uncertainty bands and the lower panels show the ratio to the results obtained for the µ 0 = H T /2 scale setting. For the transverse momentum of the hardest b-jet (p T, b 1 ), we find that the results for the three scale choices are comparable and differ in the tails only up to 5%. While the scale uncertainties for the two dynamical scales are rather constant and stay below 10%, the fixed scale setting leads to significantly larger effects of 25% or more. For the observable H had T , M b 1 b 2 and M e + µ − we find that the two dynamical scales lead to similar results. For µ 0 = µ fix larger deviations are present in the tails when comparing to µ 0 = µ dyn and µ 0 = H T /2. Specifically, the differences are up to 20% for H had T , 5% for M b 1 b 2 and 25% for M e + µ − . Furthermore, for µ 0 = µ fix the scale uncertainties increase up to even 70%-80%. For µ 0 = µ dyn and µ 0 = H T /2 they are below 20% and 10% respectively. Thus, we conclude that the two dynamical scales are alike for µ 0 . The fixed scale setting, however, is quite different already for the central value of the scale. Furthermore, it results in significantly larger scale uncertainties in the tails of these dimensionful distributions.
In figure 7 we show the azimuthal angle in the transverse plane between the two charged leptons (∆φ e + µ − ), the cosine of the angle between the two b-jets (cos (θ b 1 b 2 )), the distance in the azimuthal angle rapidity plane between the two b-jets (∆R b 1 b 2 ) and the rapidity of the Higgs boson (y H ). For ∆φ e + µ − we find that the differences between the three scales are below 3% over the entire range. The scale uncertainties are very similar for µ 0 = µ fix , µ 0 = µ dyn and µ 0 = H T /2, and only for large values of ∆φ e + µ − ≈ π, we observe a reduction of the scale uncertainties for the dynamic scale settings. In these regions, they are reduced from 18% for µ 0 = µ fix to 13% for µ 0 = µ dyn and to 9% for reduced by a special scale choice. In particular, they are reduced from 10% for µ 0 = µ fix to 7% for µ 0 = µ dyn and down to 5% for µ 0 = H T /2. This kinematical region corresponds to the case where the two b-jets (the two charged leptons) can be found almost back-to-back, which is typical for a produced top-quark pair and, to large extent, for its decay products. For ∆R b 1 b 2 and y H , we again observe small differences among the three scale settings. We notice also in this case that for ∆R b 1 b 2 , the scale dependence around ∆R b 1 b 2 ≈ π is improved for µ 0 = µ dyn and µ 0 = H T /2. For y H , the scale uncertainties are below 5% in the central rapidity region while they increase up to 15%-20% in the forward-backward regions. This observation is independent of the choice made for µ R and µ F . In conclusion, for angular distributions, which we have examined, the results are not affected by the scale choice. This is in contrast to the case of dimensionful observables, which we have discussed before. On the other hand, we find improvements in the scale dependence for dimensionless observables when we use a dynamical scale setting. In particular, for the µ R = µ F = µ 0 = H T /2 scale choice, the smallest theoretical error is obtained.
We continue the discussion at the differential level with the differential K-factor defined as K = dσ NLO /dσ LO to assess the size of the NLO QCD corrections and the perturbative stability of the two scales µ 0 = H T /2 and µ 0 = µ fix . For the discussion of the differences between a dynamical and a fixed scale, only µ 0 = H T /2 is used in the following, since the conclusions are rather similar for µ 0 = H T /2 and µ 0 = µ dyn . We start with figure 8 where the following two observables H T and p T, e + are shown. We present the LO and NLO predictions in the upper panels of the plot for µ 0 = H T /2 (left) and µ 0 = µ fix (right). The JHEP02(2022)196 lower panels display the differential K-factor and the corresponding uncertainty bands due to the scale variation. For H T , we notice that, already at LO, the scale uncertainties are reduced from 40% for µ 0 = µ fix to 32% for µ 0 = H T /2. At NLO, we find that, for the fixed scale choice, scale uncertainties can be as high as 100% and even negative values are included in the uncertainty bands. On the other hand, for µ 0 = H T /2, they are only of the order of 8%. In addition, we find that the NLO and LO predictions do not fit within the uncertainty bands in the tails for µ 0 = µ fix while for µ 0 = H T /2, the NLO predictions are completely included in the uncertainty bands of the LO prediction. The higher order QCD corrections are not constant in both cases and reach up to 35% even for the dynamical scale setting. Thus, the NLO QCD corrections are necessary for a precise prediction for H T . For the second observable, p T, e + , the dynamical scale leads to a better perturbative JHEP02(2022)196 stability in the tails where the NLO and LO predictions for µ 0 = µ fix are clearly separated from each other. The scale uncertainties are reduced from about 30% at LO to 5% at NLO for µ 0 = H T /2. Instead, for the fixed scale choice they can exceed 50% at NLO, which is even more than the corresponding LO scale uncertainty for that scale setting. Finally, the NLO QCD corrections for p T, e + with µ 0 = H T /2 are rather constant and of the order of 20% over the whole plotted range. The behavior of most of the dimensionful observables, that we have examined, is similar to that of H T and p T, e + .
In the following, we present two observables that are rather special. Specifically, in figure 9, the observables p T, b 1 b 2 and p T, miss are shown. We can notice that up to around 150 GeV for both observables, NLO QCD corrections are of the order of 20%-30%. Above this value, however, they increase substantially and can reach 90% for p T, b 1 b 2 and 60%

JHEP02(2022)196
for p T, miss . Because of these huge QCD corrections, the LO and NLO predictions do not overlap anymore for the transverse momentum values above 150 GeV. Such effects are well known and discussed (also for p T, , that is not shown here), e.g. in refs. [33,39,57], where NLO QCD calculations for top quark pair production with and without heavy bosons are presented. In ref. [92], similar findings are described for the tt production process in the di-lepton top quark decay channel at next-to-next-to-leading order (NNLO) in QCD. In this case the full NWA with corrections to top-quark production and subsequent top-quark decay is employed. Furthermore a comparison to CMS and ATLAS data is carried out. We note that such huge effects are not observed for the fix scale setting. However, since at NLO in QCD the differences between the two scale choices, µ 0 = µ fix and µ 0 = H T /2, are about 2% for p T, miss and up to 10% in the tails for p T, b 1 b 2 we can attribute these discrepancies to the corresponding LO differential cross sections. Indeed, the LO cross section, which in general is much more sensitive to the variation of the scale, changes more rapidly than the NLO curve, which drives the NLO/LO ratio in the phase-space region above 150 GeV. Observables that are constructed from the decay products of both top quarks exhibit a strong suppression of the tt cross section at LO above 150 GeV. This can be understood by analyzing tt production in the NWA. In that case the top quark decay products are boosted via their parents that have opposite transverse momenta. Thus, although the transverse momentum of the individual decay products can be quite large, the p T of the reconstructed bb, e + µ − and ν eνµ system only acquires very small values. Consequently, the transverse momenta of these systems are suppressed for p T ≥ 150 GeV. At NLO, due to additional radiation, this kinematical constraint is partially lifted resulting in an enhancement of the cross section and a large K-factor, see e.g. refs. [39,92]. For ttH production, acquiring a non-vanishing transverse momentum by the Higgs boson relaxes the kinematical constraint already at LO, which leads to much smaller but still substantial NLO QCD corrections for p T, b 1 b 2 , p T, e + µ − and p T, miss = | p T, νe + p T,νµ |, see e.g. ref. [33]. On the contrary, the scale uncertainties for p T, b 1 b 2 (p T, miss ) only decrease from 17% (12%) for µ 0 = H T /2 to 11% (10%) for µ 0 = µ fix .
At last, we show in figure 10 two examples for dimensionless observables. In particular, we display the azimuthal angle between the two b-jets, ∆φ b 1 b 2 , and the cosine of the angle between the positron and the muon, cos θ e + µ − . At LO, for both scale choices, the scale uncertainties are of the order of 30%. We find that at the NLO level in QCD for large angles, µ 0 = H T /2 leads to smaller scale uncertainties (7%) than µ 0 = µ fix (14%). For both scales, the LO and NLO predictions fit within the uncertainty bands and have similarly large QCD corrections for small angles of about 30%. For larger angles, the NLO QCD corrections are reduced to about 15% for µ 0 = H T /2 and 5% for µ 0 = µ fix . Therefore, the dynamical scale leads to smaller shape distortions, i.e. flatter differential K-factors. We conclude that, for angular distributions, both scale settings can be safely used. However, the dynamical scale choice leads to a reduction of scale uncertainties for the back-to-back phase-space region.
At the integrated level, the scale dependence clearly dominates the total theoretical uncertainties. Unquestionably, this is true at the differential level. Nevertheless, it is important to investigate the size of PDF uncertainties, since it cannot be excluded that . Differential distributions at LO and NLO QCD with the corresponding uncertainty bands for the observables ∆φ b1b2 and cos θ e + µ − for the pp → e + ν e µ −ν µ bb H process at the LHC with √ s = 13 TeV. Results are presented for µ 0 = H T /2 and µ 0 = µ fix . The NNPDF3.1 PDF set is employed. Also given are Monte Carlo integration errors. The lower panels display the differential K-factor. they can be enhanced in certain phase-space regions and play a larger role than at the integrated level. To this end, in figure 11, differential distributions for the observables H vis T , are shown for the dynamical scale setting µ 0 = H T /2 and for the following three PDF sets NNPDF3.1, CT18 and MSHT20. Each figure comprises three parts; the upper panels show the NLO prediction for three different PDF sets for µ 0 = H T /2, the middle panels display the ratio to the NNPDF3.1 set and its scale dependence, whereas the lower panels give the internal PDF uncertainties obtained for each PDF set separately. For the observable H vis T , that is defined as H vis T = p T, e + + p T, µ − + p T, b 1 + p T, b 2 , we find that the PDF uncertainties are growing quite fast towards p T tails. They are at the 5% level for CT18, 3% for MSHT20 and of the order of 2% for NNPDF3.1. For comparison, the scale uncertainties in this region are about 5%-6%. Thus, they are comparable in size

JHEP02(2022)196
to the PDF uncertainties of the CT18 PDF set. The differences between the three PDF sets are similar to the individual PDF uncertainties. The largest variation, around 3%, can be found when comparing CT18 and NNPDF3.1. For M b 1 b 2 the situation is rather similar. The PDF uncertainties increase towards the tails where they become comparable to the scale uncertainties. The latter are of the order of 8%. Also in this case, the largest differences are found between the CT18 and NNPDF3.1 PDF sets. For ∆φ b 1 b 2 the differences between the PDF sets and the PDF uncertainties are constant over the entire range. For cos (θ b 1 b 2 ) the same is true for the PDF uncertainties. However, the differences between the three PDF sets increase at larger angles. Nevertheless, theoretical predictions for all PDF sets agree well within the corresponding PDF uncertainties, that are of the order of 3% for CT18, 2% for MSHT20 and 1% for NNPDF3.1. Furthermore, they are smaller than the scale uncertainties, which are at the level of 5% for cos (θ b 1 b 2 ) and 7% for φ b 1 b 2 .
To conclude this part of the paper, the PDF uncertainties are more important at the differential than at the integrated level. They play an essential role in the tails of dimensionful observables, where they become similar in size to the corresponding scale uncertainties. For angular distributions we find that the PDF and scale uncertainties are similar in size to those at the integrated fiducial level.

Off-shell vs. on-shell top quark modelling
In this section we examine the size of full off-shell effects in the modelling of the top quarks at the integrated and differential fiducial level. So far, we have concentrated on the pp → e + ν e µ −ν µ bb H + X process at NLO in QCD. Specifically, our calculation is based on the full matrix elements for the final state e + ν e µ −ν µ bb H, where all Feynman diagrams of O(α 3 s α 5 ) are incorporated and no approximations are used. Thus, all factorizable and non-factorizable NLO QCD corrections to a 2 → 7 process are incorporated together with all the spin correlations. We note here, that non-factorizable NLO QCD corrections imply a cross-talk between production and decays of top quarks. Although we discuss Higgs production in association with a top quark pair, we take into account double-, singleand non-resonant Feynman diagrams, interferences, and off-shell effects of the top quark and the W gauge boson. The latter effects are included via incorporating Breit-Wigner propagators for the unstable particles. Consequently, the full off-shell calculation is free of ambiguities related to disentangling single-and double-resonant top-quark contributions. Such calculations, complete from the point of view of fixed-order perturbation theory, can be very time consuming. Therefore, different approximations are used to obtain results in a reasonable amount of time. A commonly used approximation to calculate processes with unstable particles is the narrow-width approximation. In this approach, the propagator of the unstable particle is expanded in the limit Γ/m → 0 leading to  TeV. The upper panels show the NLO prediction for three different PDF sets, NNPDF31, MSHT20 and CT18 and with µ 0 = H T /2. Also given are Monte Carlo integration errors. The middle panels display the ratio to the NNPDF3.1 set and its scale dependence. The lower panels give the PDF uncertainties obtained for each PDF set separately.

JHEP02(2022)196
into on-shell production and on-shell decay. This approximation reduces the complexity of the process enormously, as only double resonant diagrams are taken into account, whereas all single-and non-resonant diagrams and their interferences are systematically neglected. In our case, we treat the two top quarks and the two W -bosons in the NWA, which leads to the following decay chain In this case, the Higgs boson can only be produced in the production of the top quark pair, because m t < m W +m H . With both calculations at our disposal, we can examine the range of applications of theoretical predictions in NWA at the integrated and differential fiducial level. In addition to the NWA and the full off-shell calculation, we will also present the results for a NWA, in which no QCD corrections in the decay of the top quarks are included. We label the latter NWA LOdec . By comparing the NWA and NWA LO dec predictions, on the other hand, we can examine how the inclusion of QCD corrections in the top-quark decays affects our fixed-order NLO calculation at O(α 3 s α 5 ). 5 In the NWA case, the top quark width as given in eq. (3.5) is used. Specifically, we employ Γ NLO t, NWA for the NWA results and Γ LO t, NWA for NWA LOdec at NLO. Furthermore, we do not expand the top quark width in a series in the strong coupling constant in the NWA calculation, since this is not possible in the full off-shell calculation and therefore would affect our study of the size of the non-factorisable corrections for the process. The calculations for the NWA and NWA LOdec are also performed within the Helac-Nlo framework, that has recently been extended to provide theoretical predictions in this approximation [66]. The Catani-Seymour dipole formalism is used, and the phase-space restriction on the dipoles is employed to check the independence of the real emission results on this parameter.
In table 4, the integrated cross section at LO and NLO QCD is shown for µ 0 = H T /2 and µ 0 = µ fix . Results are obtained with the default PDF set NNPDF3.1. The full off-shell calculation, NWA and NWA LOdec results are presented with the corresponding Monte Carlo error and a theoretical uncertainty estimate obtained through scale variation. We note here, that in the full off-shell calculation the scale uncertainties are calculated with the 7-point scale variation while for NWA and NWA LOdec the scale uncertainties are obtained with the 3-point scale variation. The latter approach, however, leads to similar scale uncertainties as the former one, as discussed in section 5. At LO, the deviations between the NWA and the full off-shell calculation are about 0.5% for µ 0 = H T /2 and 0.3% for µ 0 = µ fix . Thus, they are negligible compared to the large scale uncertainties of about 30%. At NLO, the off-shell effects are about 0.4% for both scale settings and, again, significantly smaller than the scale uncertainties of about 5%. These differences are consistent with the expected error of the NWA that is of the order of Γ t /m t ≈ 0.8% [93]. We notice that the maximum of the scale variation for both scales is essentially reproduced in the NWA and that the differences between the full off-shell calculation and the NWA do not exceed 0.5% for all scale choices required in the 3-point scale variation.  in the NWA LOdec is about 5% (6%) larger for µ 0 = H T /2 (µ 0 = µ fix ) than the result in the NWA. Thus, the inclusion of QCD corrections in the top-quark decays leads to a reduction of the integrated cross section. Moreover, the scale uncertainties increase to about 9% when the NLO QCD corrections in the top-quark decays are omitted.
Even though full off-shell effects of the top quarks and W gauge bosons are formally suppressed by the corresponding widths and are small for the inclusive cross section, they can be strongly enhanced in exclusive observables. Specifically, for dimensionful observables at the differential level, the contribution from single-and non-resonant Feynman diagrams can be enhanced and even become larger than the contribution from double resonant diagrams in certain phase-space regions. The latter comprise high p T regions and kinematical edges, that are related to the masses of the on-shell particles and the cuts that all final states have to pass. As an example, in figure 12, we show differential distribution at NLO in QCD for the observables M (be + ) min , H T , M b 1 b 2 and p T,b 1 b 2 . We provide theoretical results with the full off-shell effects included as well as predictions in the NWA and NWA LOdec . We utilize µ 0 = H T /2 and the NNPDF3.1 PDF set. The lower panels display the ratios of the corresponding NWA result to the full off-shell one. We additionally provide theoretical uncertainties as obtained from the scale dependence for the full off-shell case. The observable M (be + ) min is defined as the minimum mass of the positron and a b-jet. At LO in the NWA, this observable is limited by m 2 t − m 2 W ≈ 153 GeV due to the finite mass of the top quarks and W -bosons. At NLO, this kinematical edge is smeared out by real radiation, but the NWA is still unable to recover the full calculation in the vicinity of M (be + ) min ≈ 153 GeV. For M (be + ) min ≤ 153 GeV, the NWA leads to good agreement with the off-shell calculation and the differences between the two cases are significantly smaller than the scale uncertainties. The missing corrections in the decay of the top quark lead to further shape distortions in these regions of the phase space, which, this time however, are larger than the scale variation of the full result. For M (be + ) min ≥ 153 GeV, on  the other hand, we obtain substantial non-factorizable corrections that are of the order of 50%-60%. For H T , we observe that the differences between the NWA and the off-shell calculation are much smaller than the scale uncertainties at the beginning of the spectrum. Towards the tails, these differences rise up to about 10%, which is roughly the size of the corresponding scale uncertainties. In the tails, NWA and NWA LOdec are very similar, while for H T ≈ 500 GeV we find additional shape distortions for the latter approach, close to 10% . For M b 1 b 2 and p T, b 1 b 2 we once more observe that the NWA and the off-shell calculation agree well for small and moderate values. However, in the tails, the differences are equal to the scale uncertainties of about 10% for M b 1 b 2 and 20% for p T, b 1 b 2 . Also here the NWA LOdec leads to additional shape distortions of about 10%.  In figure 13, we exhibit differential distributions for the two dimensionless observables ∆φ b 1 b 2 and y e + . For ∆φ b 1 b 2 , we find that the NWA is able to recover the full off-shell calculation over the entire range. Thus, full off-shell effects are negligible compared to scale uncertainties. On the other hand, NWA LOdec leads to shape distortions in the range of 2%-10%. The latter cannot be simply corrected by a re-scaling factor. Likewise, for the rapidity of the positron, y e + , we find that full off-shell effects are rather insignificant. Furthermore, the NWA LOdec approach introduces an overall normalisation correction of the order of +5%.
Concluding this part, full off-shell effects can be safely neglected at the integrated level, while they might be substantial for precise predictions at the differential level. For dimensionless observables, the off-shell effects are negligible, but for dimensionful observables close to kinematical edges or in the high-energy tails of various dimensionful distributions, these effects become significant. They are either similar in size or well in excess of the corresponding scale uncertainties. Consequently, they cannot simply be omitted. When studying differential cross sections for the process at hand, a dedicated and careful analysis of such effects is required. We further add that the NWA LOdec does not provide the correct normalisation at the integrated fiducial level. At the differential level, already for some angular distributions, additional shape distortions from the missing NLO QCD corrections in the top quark decays are observed. We also find additional effects for many dimensionful observables over the entire range. These shape distortions are usually larger than the scale uncertainties of the full off-shell calculation.

JHEP02(2022)196 8 Bottom quarks in the initial state
The last point in the discussion of the pp → e + ν e µ −ν µ bb H process with a stable Higgs boson that we would like to scrutinise concerns the contribution of initial state bottom quarks. Such contributions are often neglected in similar calculations due to the suppressed bottom quark PDFs. We investigate the size of these contributions in comparison to the other theoretical uncertainties discussed in the sections 5 and 6. The inclusion of bottom quarks in the initial leads at NLO to additional sub-processes such as bg → e + ν e µ −ν µ bbb H andbg → e + ν e µ −ν µ bbb H, in which we have three bottom quarks in the final state. As we are working in the five-flavour scheme were the b quarks are treated as massless, such sub-processes have an additional infrared singularity because of the gluon splitting into a bottom quark pair, which has to be resolved. This is possible if we take into account the flavour of the bottom quarks in the jet algorithm and require that two b-jets with opposite charge are recombined into a gluon bb → g, such that these events do not pass the event selection. Such collinear divergences do not affect the recombination of two b partons with the same charge. Thus, for bb andbb we can choose which definition we would like to apply. In particular, we use two different variants as described in ref. [94].
In the first scheme, that we label charge-blind b-jet tagging, the charge of the b-jets is neglected. We treat two b partons with the same charge in the same way as a bb pair and also recombine them into a gluon. Therefore, we have the following recombination rules for this scheme: Since the charge is neglected, we require that at least two b-jets with arbitrary charge pass the cuts outlined in section 3, otherwise the event is rejected. Consequently, already at LO, besides the sub-process bb → e + ν e µ −ν µ bb H, we must also include the following two additional sub-processes bb → e + ν e µ −ν µ bb H andbb → e + ν e µ −ν µbb H. (8.2) In the second scheme, that we label charge-aware b-jet tagging, we no longer neglect the charge of the b-jet constituents and require the following recombination rules We further require that at least two b-jets with opposite charge pass the cuts. The two new sub-processes, which are necessary in the charge-blind scheme, do not pass this requirement and can be safely neglected. From an experimental point of view, the additional charge tagging of the b-jets might lead to smaller b-tagging efficiency and larger statistical uncertainties.
In order to validate our results for both schemes, the computation has been performed for two extreme cases of the phase-space restriction parameter of the subtraction terms in the Nagy-Soper subtraction scheme. In this way, we check the independence of the real emission part on this parameter. For each case, a perfect agreement has been found. 0.26% 0.33% Table 5. Integrated fiducial cross section at LO and NLO QCD for the pp → e + ν e µ −ν µ bb H process at the LHC with √ s = 13 TeV. Results are shown for µ 0 = H T /2 and µ 0 = µ fix . The NNPDF3.1 PDF set is employed. Predictions are given without and with initial state bottom quark contributions in the charge-aware and charge-blind scheme. Last two columns present the ratio δ X = (σ X − σ no b )/σ X .

JHEP02(2022)196
In table 5, the integrated fiducial cross sections at LO and NLO QCD are given without (σ no b ) and with the initial state bottom quark contribution in the charge-aware and chargeblind scheme (denoted as σ aware and σ blind respectively). Also shown is the ratio δ X = (σ X − σ no b )/σ X , where X = aware , blind. Furthermore, Monte Carlo errors and scale uncertainties are reported for each case. We use the two scale settings µ 0 = H T /2 and µ 0 = µ fix , and our default PDF set NNPDF3.1. As mentioned before, our dynamical scale setting does not require top-quark reconstruction, and, therefore, a precise estimate of the bottom-initiated contribution can be straightforwardly obtained. On the other hand, for µ 0 = µ dyn , an additional procedure would have to be introduced to correctly assign the two b-jets to the corresponding top quarks. At LO, we find that the initial state bottom quark contribution amounts to about 0.2% independently of the scale setting and the scheme employed. Indeed, the contribution of the additional bb andbb channels in the charge-blind scheme is below the integration error of about 0.01%. Compared to the scale uncertainties of about 30%, the bottom-initiated contribution can be completely neglected at LO for the process at hand. At NLO, we find that the initial state bottom quark contribution increases due to the additional gb/gb channels with up to three b-jets in the final state. However, the overall impact is still very small. Specifically, in the charge-blind scheme the bottom-initiated contribution accounts to 0.3% for both µ 0 = H T /2 and µ 0 = µ fix . In the charge-aware scheme we have 0.2% for µ 0 = H T /2 and 0.3% for µ 0 = µ fix . The bottominitiated contribution that is more than fifteen times smaller than the corresponding scale uncertainties, has no phenomenological impact and can be safely omitted. As this is a negligible fraction of the integrated fiducial cross section, the choice of a particular scheme is irrelevant. On the other hand, when only sub-processes with initial state bottom quarks are considered, theoretical predictions for the charge-blind scheme are about 25%-30% larger than the corresponding ones for the charge-aware scheme.
We also conduct a similar comparison at the differential level to check whether the bottom-initiated contribution is enhanced in some particular regions of the phase space. We expect that the largest impact, if any, should be found in hadronic observables due to   2 and cos (θ b 1 b 2 ). The transverse momentum of the Higgs boson is also given in figure 14 because we would like to study the impact of the bottom-initiated sub-processes on this important observable. All results are obtained for µ 0 = H T /2 and the NNPDF3.1 PDF set. Indeed, we observe that for p T,b 1 and M b 1 b 2 the initial state bottom quark contribution increases towards the tails relative to the calculation without this contribution. Specifically, we find effects up to about 3%. Compared to the asymmetric scale uncertainties the initial state bottom quark contribution is similar in size, although, still smaller than the conservative scale uncertainty estimate of about 10%. The latter can be obtained by selecting the maximum of the two results in the bin that represents the theoretical error for the observable under consideration. Moreover, we find no significant differences between the two tagging schemes.
For hadronic angular distributions such as cos (θ b 1 b 2 ), the differences between the calculations do not exceed 0.6% at large angles and decrease even further towards small angles. For large angles, we again find that the initial state bottom quark contribution is similar in size to the asymmetric scale uncertainties but significantly smaller than the conservative scale uncertainties of about 5%. Finally, for non-hadronic observables such as the transverse momentum of the Higgs boson we find no significant effects.
To sum up this part of the paper, the bottom-initiated contribution is negligible at the integrated and differential level when comparing to theoretical uncertainties as estimated by varying the renormalization and factorization scales in α s and PDFs. However, due to the fact that in the tails of some dimensionful hadronic observables these contributions can increase up to 3%, they can be similar in size or even larger than the corresponding internal PDF uncertainties. For dimensionless hadronic and non-hadronic observables the initial state bottom quark contribution does not play an important role.

Theoretical predictions with Higgs boson decays
In the last part of the paper we discuss the cross section at the integrated and differential level but this time various decay channels of the Higgs boson are included as well. As it is well known, the Higgs boson can only be detected via its decay products. Thus, it is important to obtain precise theoretical results for fully realistic final states. In our framework, this can be done without further time consuming calculations, because we simply reuse the LHEFs generated for the pp → e + ν e µ −ν µ bb H process and model the decay of the Higgs boson in the NWA afterwards.
In table 6, integrated cross sections are provided for the pp → e + ν e µ −ν µ bb H process with Higgs boson decays into bb, τ + τ − , γγ and e + e − e + e − . For comparison, predictions JHEP02(2022)196 with a stable Higgs boson have also been included. All results are obtained for µ 0 = H T /2 with our default NNPDF3.1 PDF set. For the decay into a bottom quark pair, in which NLO QCD corrections are taken into account, we also use the following renormalisation scale setting µ R, dec = m H in the decay. In addition, Monte Carlo errors and scale uncertainties are presented. Also given in the last column is the K factor. For all cases, except for the decay into a bottom quark pair, the 7-point scale variation is used. For H → bb a 21-scale variation is used instead, where µ R, dec is varied independently of the renormalisation and factorisation scales. The cross sections of the different decay channels are smaller than the naive multiplication of the production cross section with the branching ratio. This is of course due to the additional cuts that are applied on the Higgs boson decay products. We can estimate the number of events in the LHEFs that pass the event selection by the ratio of these two values. For all 1 → 2 decays approximately 70%-80% of events pass the event selection, while for the 1 → 4 decay chain only about 15% of events remain. Thus, the cuts on the final state of the last decay mode have a greater effect on the kinematics of this process and also lead to larger statistical uncertainties due to the extremely limited available phase space. For Higgs decays into τ + τ − and γγ, we find that the integrated cross section behaves very similarly to the case with a stable Higgs boson. The scale uncertainties at LO and NLO and the K-factor are very similar, since both decay channels differ only slightly in the overall kinematics due to the cuts. In contrast, for the Higgs boson decay into e + e − e + e − , we observe that the theoretical uncertainties increase from about 5% to 6% and that the K-factor rises from 1.23 to 1.30. The differences between these three decay channels, however, are rather small compared to the Higgs boson decay into a bb pair. For the latter case, we find already at LO that the scale uncertainties increase from 30% to 44%. The growth is caused by the running of the Yukawa coupling in the MS scheme that is taken into account in the scale variation. At NLO, the scale uncertainties drop from 44% to about 10%, but are still a factor of 2 larger than the corresponding scale uncertainties as estimated for the case of a stable Higgs boson. In addition, for the H → bb decay, the K-factor decreases from 1.23 to 1.14 leaving us with smaller QCD corrections.
To investigate the size of NLO QCD corrections in the Higgs boson decay, we again present the integrated fiducial cross section at NLO for the pp → e + ν e µ −ν µ bb H process, but this time the H → bb decay is included only at LO. Our NLO QCD result, as obtained with the NLO Yukawa coupling, is

JHEP02(2022)196
At the differential level, we first discuss the pp → e + ν e µ −ν µ bb H process with a subsequent H → bb decay. Since the production and decay of the Higgs boson are disentangled, the QCD corrections in the Higgs boson decay do not affect substantially non-hadronic observables except for the overall normalisation. On the other hand, we expect additional shape distortions for hadronic observables that can be constructed from the four b-jets. Indeed, the latter can now originate from the production or from the decay of the Higgs boson. For example, in the case of b-jets originating from top quark decays, these shape distortions can be caused by the possible recombination of the bottom quark from t → bW + ort →bW − decay with a gluon emitted in H → bb. Furthermore, an additional gluon in the Higgs boson decay, if resolved as an extra light jet by the jet algorithm and independently of whether it is observed or not, can affect the reconstruction of the Higgs boson system. Such effects can impact various b-jet distributions mainly in the soft and collinear regions of the phase space. Thus, we might expect some modifications both in dimensionful and dimensionless observables.
In figure 15, we show differential distributions for the following observables p T, b 1 , for the pp → e + ν e µ −ν µ bb H (H → bb) process at the LHC with √ s = 13 TeV. We use our default NNPDF3.1 PDF set and employ the following scale settings µ 0 = H T /2 and µ R, dec = m H . Results are again presented at LO, NLO and for the NLO LOdec H case. The lower panels display the ratio to the latter case. We distinguish between b-jets coming from the top quarks (b 1 , b 2 ) and the Higgs boson (b H 1 , b H 2 ) via the superscript H. The Higgs boson reconstruction is performed as described in section 3. For the transverse momentum of the hardest b-jet coming from top quark decays, denoted as p T,b 1 , we find that the scale uncertainties in the tails are substantially reduced when NLO QCD corrections are included. Indeed, they decrease from about 50% at LO to 20% at NLO LOdec H . The inclusion of higher-order corrections in the decays of the Higgs boson additionally influences these theoretical uncertainties, reducing them even further to about 15%. All three results are in good agreement within the corresponding uncertainty bands.
For small values of p T,b 1 an enhancement in the range of 30%-50% can be noticed. The latter is caused by gluon emission in the Higgs boson decay that affects the Higgs boson reconstruction. On the other hand, in the tails of this distribution, we find only small corrections of the order of 10%. For the invariant mass of two b-jets originating from the tt pair, we observe a very similar reduction in scale uncertainties when going from LO to NLO. Again, all results agree well within the corresponding scale uncertainties. Moreover, at the beginning of the spectrum, the inclusion of the NLO QCD corrections in the Higgs boson decay increases the cross section by about 25%, while we observe no significant corrections in the tails. Looking at the angular separation between the two b-jets, denoted as ∆R b 1 b 2 , we find that, close to ∆R b 1 b 2 ≈ 3, the scale uncertainties decrease from 45% at LO to 15% at NLO LOdec H and 10% at NLO. This phase-space region corresponds to the characteristic back-to-back production of the two b-jets from the top-quark pair. On the other hand, for small values of ∆R b 1 b 2 , we find an enhancement up to even 30%, which clearly shows that these corrections are caused by the mismatch in the Higgs boson reconstruction. We plan to investigate this issue in more detail in the near future. To this end, different reconstruction techniques, that mimic experimental analyses as closely as possible, will be employed and their impact on various observables will be investigated. Finally, in the case of the cosine of the two b-jets originating from the Higgs boson decays, denoted as cos , the scale uncertainties decrease from about 45% at LO to 15%-20% at NLO LOdec H . The additional NLO QCD corrections in the Higgs boson decay lead to a further reduction of theoretical uncertainties down to about 10%. Furthermore, these corrections are below 5% at large angles and of the order of 10% for small values of In a next step, we discuss the differences between the various Higgs boson decay channels at the differential level. For this purpose, we examine the normalised differential cross-section distributions to assess whether only the overall normalisation is different, or whether we can find further shape distortions. In figure 16 we present results for the 1 d /dp T, b1 [ Figure 16. Normalised differential distributions for the observables p T, b1 , |∆y e + µ − |, p T,H and cos (θ b1b2 ) for the pp → e + ν e µ −ν µ bb H process at the LHC with √ s = 13 TeV. Results are given for µ 0 = H T /2 and µ R, dec = m H . The NNPDF3.1 PDF set is employed. The case of the stable Higgs is displayed as well as predictions with the Higgs boson decays into bb, τ + τ − , γγ and e + e − e + e − . Also given are Monte Carlo integration errors. The lower panels display the ratio to the case with the stable Higgs boson.
pp → e + ν e µ −ν µ bb H process with a subsequent Higgs boson decay into bb, τ + τ − , γγ or e + e − e + e − . For comparison, also the case with a stable Higgs boson is displayed. We present the following observables p T, b 1 , |∆y e + µ − |, p T,H and cos (θ b 1 b 2 ). The lower panels display the ratio to the case with a stable Higgs boson. We use the same scale settings and PDF set as for the integrated fiducial cross section. As in the previous discussion, we JHEP02(2022)196 distinguish the decay products from the top quark and Higgs decays by using the Higgs boson reconstruction procedure outlined in section 3. We use the Higgs boson reconstruction for all decay channels, but for the decay channels without bottom quarks the MC truth and the Higgs boson reconstruction is exactly the same due to the absence of additional electroweak radiation. For the transverse momentum of the hardest b-jet, we find that the decay into four electrons leads to a completely different shape. This is not surprising, of course, since we are talking about the decay pattern 1 → 4. In particular, the p T,b 1 spectrum is significantly harder due to the strict event selection. This observation essentially holds for all dimensionful observables. By comparison, our results for all 1 → 2 decay modes are quite similar. The NLO QCD corrections in the H → bb decay lead to an enhancement at the beginning of the spectrum up to 40%. However, in the tails, the decay channel bb differs only by about 10% from the case with a stable Higgs boson. For the remaining two decay channels, τ + τ − and γγ, differences up to 2% are only observed. For the rapidity separation of the two charged leptons |∆y e + µ − |, we find that all results are very similar for small values. If the two leptons are well-separated in rapidity, the differences increase, especially for the Higgs boson decay into four electrons. For the 1 → 2 decay channels, discrepancies are below 10%. Indeed, in the case of H → bb and H → γγ we find differences up to 8%, which are reduced to 3% for the H → τ + τ − case. For the transverse momentum of the Higgs boson, we find differences up to about 10% for γγ, bb and for τ + τ − . Thus, the inclusion of QCD corrections in the Higgs boson decay into bb does not lead to significant enhancements in any specific phase-space region. However, these differences are comparable in size to the scale uncertainties and thus can not be simply ignored. We note, of course, that for the H → e + e − e + e − decay chain, the p T,H spectrum has a completely different shape. Last but not least, for cos (θ b 1 b 2 ) we observe for bb once again that the region of small angles is enhanced by up to 15% due to the mismatch in the reconstruction of the Higgs boson. Similar effects in absolute value can be observed for the H → e + e − e + e − decay chain. We find that the other two decay channels lead to similar results, which still differ by about 7%-8% from the prediction for a stable Higgs boson. Also in this case, these shape distortions are similar in size to the scale uncertainties. For large values of θ b 1 b 2 , smaller discrepancies well below 10%, can be noticed for H → bb and H → e + e − e + e − . They are again similar in absolute magnitude but opposite in sign.

Summary
In this paper, we discussed the production and decay of the Higgs boson in association with a tt pair. In particular, we calculated NLO QCD corrections to pp → e + ν e µ −ν µ bb H (H → X) at the LHC with √ s = 13 TeV including full off-shell effects of the top quarks and W gauge bosons. In the computation these unstable particles have been described by Breit-Wigner distributions. Furthermore double-, single-and non-resonant top quark/W gauge boson contributions as well as all interference effects have been consistently incorporated already at the matrix element level. The calculation of the production process was performed within the Helac-Nlo framework. The decays of the Higgs boson were modelled in the NWA with the help of the LHEFs.

JHEP02(2022)196
For Higgs production, we first compared our results with the findings from ref. [33] both at the integrated and differential level. For the two scale settings, µ 0 = µ dyn and µ 0 = µ fix , that were employed there, we found an overall good agreement between the two calculations. In the next step, we introduced an additional scale choice, µ 0 = H T /2, and performed a comparison between the three scales. At the integrated fiducial level, we found that, for all scale settings, the NLO QCD corrections are of the order of 20%. The scale uncertainties decreased from about 30% at LO to 5% at NLO. All results are in good agreement within the corresponding scale uncertainties. Furthermore, we reported the PDF uncertainties for the following PDF sets NNPDF3.1, CT18 and MSHT20. These PDF uncertainties are up to only 1%-2%. Thus, the theoretical uncertainties due to the scale dependence are the dominant source of the theoretical systematics.
A similar comparison was performed at the differential level. We found that the fixed scale setting led to perturbative instabilities for dimensionful observables in high-energy tails. Even at the level of differential cross sections, the fixed scale choice resulted in larger scale uncertainties. The differences between the other two dynamical scales are only minor. For µ 0 = H T /2 we obtained NLO QCD corrections up to 20%-30% and scale uncertainties of the order of 5%-10%. In addition, the PDF uncertainties can increase in the tails of dimensionful observables. Specifically, for the CT18 PDF set we observed PDF uncertainties of about 5%-7%. The latter are similar in size to the corresponding scale uncertainties. However, for our default NNPDF3.1 PDF set, the PDF uncertainties decreased to 2%-3% only. Therefore, for our default setup, also at the differential level, the theoretical error is dominated by the theoretical uncertainties resulting from the scale variation.
In the next step, we discussed the off-shell effects of the top quarks by comparing the full off-shell calculation with the results in the NWA and NWA LOdec at the integrated and differential level. At the integrated level, the off-shell effects are of the order of 0.3%-0.5%. The missing NLO QCD corrections in the decays of the top quarks increase the integrated fiducial cross section by about 5% and the scale uncertainties to 9%. At the differential level, the full off-shell top quark effects are significant close to kinematical edges and in high energy regions of the phase space. In the tails of dimensionful observables the full off-shell effects are in the range of 10%-20% and thus similar in size to the scale uncertainties. In the case of special observables, like for example M (be + ) min , they can be as high as 50%-60%. For angular distribution, the off-shell effects do not play a crucial role. In addition, the NWA LOdec leads to a wrong normalisation and further shape distortions in dimensionful observables and even angular distributions.
Subsequently, we studied the bottom-induced contributions to the process at hand. The latter is often neglected in similar calculations due to the suppressed bottom quark PDF. We used two different schemes for b-jet tagging, the so-called charge-blind and chargeaware tagging of b-jets, to obtain IR-safe results. At the integrated level the contribution is negligible of the order of 0.2%-0.3% only. At the differential level the contribution from the initial state bottom quarks can be enhanced in the tails of dimensionful hadronic observables, reaching up to 3%. The two different schemes lead essentially to the same phenomenological results.

JHEP02(2022)196
Last but not least, we combined ttH production with Higgs boson decays. We considered H → bb, H → τ + τ − , H → γγ and H → e + e − e + e − . At the integrated fiducial level, we observed that the scale dependence and the K-factor are very similar to the case of a stable Higgs boson for all decays except for H → bb. For this decay channel, we included NLO QCD corrections in the decay and performed the renormalisation of the bottom Yukawa coupling in the MS-scheme, which led to a scale-dependent Yukawa coupling. Consequently, the scale uncertainties are about 45% at LO and of the order of 10% at NLO. The K-factor is 1.14 and, therefore, smaller than the one for the case of a stable Higgs boson. This reduction is caused by the smaller Yukawa coupling of the bottom quark at NLO. Without corrections in the decay the scale uncertainties increased to 14% and the cross section decreased by about 5%. At the differential level the scale dependence is reduced by about 5% due to the NLO QCD corrections in the decay. We obtain large shape distortions of about 25% for b-jet observables at lower values of dimensionful observables and for small opening angles due to the Higgs reconstruction. On the other hand, the Higgs boson decay into four electrons leads to completely different normalised distributions. The latter effect can be especially visible for dimensionful observables, where significantly harder spectra are observed due to the strict cuts on the four charged leptons. For the Higgs boson decays into τ + τ − and γγ in normalised differential distributions shape distortions up to 14% have been found, compared to the case with a stable Higgs boson.
To recapitulate, non-factorisable NLO QCD corrections as well as higher-order QCD effects in top-quark decays impacted significantly the pp → e + ν e µ −ν µ bb H cross section in various phase-space regions. Furthermore, realistic predictions for different Higgs boson decay modes showed significant distortions compared to the case with a stable Higgs boson. Therefore, to obtain accurate predictions, all these effects should be taken into account in future comparisons between theoretical results and experimental data.