NLO QCD corrections to off-shell top-antitop production with leptonic decays in association with a Higgs boson at the LHC

We compute the hadronic production of top-antitop pairs in association with a Higgs boson at next-to-leading-order QCD, including the decay of the top and antitop quark into bottom quarks and leptons. Our computation is based on full leading and next-to-leading-order matrix elements for $e^+ \nu_e \mu^-\bar{\nu}_\mu b \bar{b} H(j)$ and includes all non-resonant contributions, off-shell effects and interferences. Numerical results for the integrated cross section and several differential distributions are given for the LHC operating at 13 TeV using a fixed and a dynamical factorization and renormalization scale. The use of the dynamical instead of the fixed scale improves the perturbative stability in high-energy tails of most distributions, while the integrated cross section is hardly affected differing by only about one per cent and leading to the same K factor of 1.17.


Introduction
After the discovery of a new boson of a mass around 125 GeV by the CMS and ATLAS Collaborations [1,2] at the LHC, the determination of its quantum numbers and couplings to other particles has become a high priority in particle physics. Results from the first run of the LHC strongly support the hypothesis that this particle is the Higgs boson predicted by the Standard Model (SM) of particle physics. In the Brout-Englert-Higgs symmetry breaking mechanism of the SM the Higgs boson is the key for understanding the origin of mass. Since in this framework the Higgs boson couples to fermions with a strength proportional to their mass via Yukawa interactions, its coupling to the heaviest quark, the top quark, is of particular interest. The main production mechanism of the Higgs boson in the SM is gluon fusion, gg→H. This process is sensitive to the top-quark Yukawa coupling, but possible heavy particles beyond the SM running in the loop bias its determination. On the other hand, the production of a SM Higgs boson in association with a top-quark pair allows a direct access of the top-quark Yukawa coupling already at tree level, disentangling it from possible beyond-SM contributions.
Several searches for the associated production of the Higgs boson with a top-quark pair for a variety of decay channels have been published by ATLAS [30][31][32][33][34] and CMS [35][36][37][38][39][40]. These searches are challenging due to a large background from ttbb and ttjj production, and so far no evidence of a ttH signal over background has been found. The current ratio of the measured ttH signal cross section to the SM expectation quoted by ATLAS is µ = 1.5 ± 1.1 [34] and by CMS µ = 1.2 +1.6 −1.5 [40] for a Higgs-boson mass of 125 GeV. In this article we present the calculation of the NLO QCD corrections to the hadronic production of a positron, a muon, missing energy, two b jets and a SM Higgs boson, which we assume to be stable (pp → e + ν e µ −ν µ bbH) at the 13 TeV LHC, which includes the resonant production of a top-antitop-quark pair in association with a Higgs boson with a subsequent leptonic decay of the top and the antitop quark. Our calculation includes all NLO QCD correction effects in ttH production and top decays and also takes into account all off-shell, non-resonant and interference effects of the top quarks. We consider the fixed renormalization and factorization scale used in Ref. [9] and alternatively the dynamical scale choice from Ref. [13] and investigate their quality in reducing the dependence of the integrated cross section as well as differential distributions on the factorization and renormalization scales. The phase-space integration is performed with a newly implemented in-house multichannel Monte Carlo program, using phase-space mappings similar to Ref. [41]. The Monte Carlo implements the dipole subtraction method [42][43][44][45] for the computation of the real corrections and is linked to the matrix-element generator Recola [46] for the computation of the LO and NLO matrix elements as well as colour-and spin-correlated squared matrix elements needed for the evaluation of subtraction terms.
The calculation follows in many respects the one of pp → e + ν e µ −ν µ bb in Ref. [47]. In particular, since the additional Higgs boson in the final state has no colour charge, the contributing Catani-Seymour dipoles [42,44] are the same. For the hadronic process pp → e + ν e µ −ν µ bbH we find a qualitative similar behaviour of NLO corrections to the integrated cross section as well as for differential distributions of the same observables. Where appropriate we compare results and point out differences. Moreover, we compare our results for the total cross sections to those of existing calculations for on-shell Higgs boson and top quarks [9,13].
The paper is organised as follows. In Section 2 we discuss details of the calculation such as contributing subprocesses and Feynman diagrams and explain some technical aspects of the real (Section 2.1) and virtual (Section 2.2) corrections. We present numerical results for the LHC operating at √ s = 13 TeV in Section 3. In Section 3.1 we list and explain the input parameters, jet definition, cuts and scale choices for our calculation, while results for integrated cross sections and a discussion of the scale dependence are provided in Section 3.2. In Section 3.4 we display and discuss several differential distributions for a fixed and a dynamical scale choice. We have performed several checks of our calculation which we present in Section 4. Our conclusions are given in Section 5.

Details of the calculation
We compute the QCD corrections to the full hadronic process pp → e + ν e µ −ν µ bbH. (2.1) We consider the tree-level amplitude at O α s α 5/2 including all resonant, non-resonant, and off-shell effects of the top quarks and all interferences. Neglecting flavour mixing as well as contributions from the suppressed bottom-quark parton densities and counting u, d, c and s quarks separately, we distinguish 5 partonic channels for the LO hadronic process: the gluon-induced process gg → e + ν e µ −ν µ bbH and four processes from qq → e + ν e µ −ν µ bbH by substituting different quark flavours (q = u, d, c, s). Throughout this paper we consider the bottom quark massless, implying no contribution from tree diagrams involving the Higgs-bottom-quark coupling. The gg process involves 236 and the qq processes 98 tree diagrams each under these prerequisites. In Figure 1 we show sample diagrams grouped by the number of top-quark resonances.  Figure 2: Representative hexagon and heptagon one-loop Feynman diagrams with two top-quark resonances.

Real corrections
The real correction process pp → e + ν e µ −ν µ bbHj receives contributions from the 13 partonic subprocesses gg → e + ν e µ −ν µ bbHg, qq → e + ν e µ −ν µ bbHg, gq → e + ν e µ −ν µ bbHq, gq → e + ν e µ −ν µ bbHq, where the gg process involves 1578 tree diagrams and the qq, gq and gq processes, all related by crossing symmetry, 614 tree diagrams each. Gluon Bremsstrahlung in the real corrections gives rise to IR divergences by soft or collinear configurations, which cancel for the final state for infrared-safe observables upon combination with the virtual corrections. Singularities from collinear initial-state splitting factorize and can be removed by MS redefinition of the parton distribution functions. We employ the Catani-Seymour subtraction formalism [42,44] for the regularization and analytical cancellation of IR singularities. Both the amplitudes for the real-correction subprocesses as well as the colour and spin-correlated amplitudes of the subtraction terms have been calculated with Recola.

Virtual corrections
The partonic subprocesses for the virtual QCD corrections can be identified with those at LO. We compute the virtual corrections in the 't Hooft-Feynman gauge, where the gg process involves 9074 loop diagrams and the qq processes 2404 loop diagrams each. The most complicated one-loop diagrams are heptagons (Sample diagrams are displayed in Figure 2).
The resonant top quarks, Z bosons and W bosons are treated in the complex-mass scheme [48][49][50], where the masses of unstable particles are consistently treated as complex quantities leading in particular to a complex weak mixing angle, For the renormalization we use the on-shell renormalization scheme as described in Ref. [49] for the complex-mass scheme.
For the computation of the matrix elements for the virtual corrections we employ Recola [46] in dimensional regularisation, which integrates the Collier [51,52] library for the numerical evaluation of one-loop scalar [53][54][55][56] and tensor integrals [57][58][59]. We compared our results for the virtual NLO contribution to the squared amplitude, 2 Re M * 0 M 1 , for many phase-space points with MadGraph5_aMC@NLO [60] (see Section 4 for details).

Input parameters, jet definition, cuts and scale choice
We present results for integrated cross sections and differential distributions for the LHC operating at √ s = 13 TeV. For the computation of the hadronic cross section we employ LHAPDF 6.05 with CT10NLO parton distributions at LO and NLO QCD. We use the value of the strong coupling constant α s as provided by LHAPDF based on a one-loop (two-loop) accuracy at LO (NLO) with N F = 5 active flavours. In the renormalization of α s the topquark loop in the gluon self-energy is subtracted at zero momentum. The running of α s in this scheme is generated by contributions from light-quark and gluon loops only. For the fixed renormalization and factorization scale µ fix = 236 GeV, defined below in (3.10), we find We neglect contributions from the suppressed bottom-quark parton density. The electromagnetic coupling α is derived from the Fermi constant in the G µ scheme [61], We compute the width of the top quark Γ t for unstable W bosons and massless bottom quarks according to Ref. [62]. At NLO QCD it reads For the top-quark width at LO we neglect the correction term −2α s F 1 (y)/(3π).
As input, we employ the following numerical values for the masses and widths: which includes the measured value of the Higgs-boson mass with zero width since we assume it to be stable. We neglect the masses and widths of all other quarks and leptons.
We convert the measured on-shell values (OS) for the masses and widths of the W and Z boson into pole values for the gauge bosons (V = W, Z) according to Ref. [63], that enter the calculation. We use the anti-k T algorithm [64] for the jet reconstruction with a jet-resolution parameter R = 0.4. The distance between two jets i and j in the rapidity-azimuthal plane is defined as with the azimuthal angle φ i and the rapidity y i = 1 2 ln E+pz E−pz of jet i, where E is the energy and p z the component of momentum along the beam axis. Only final-state quarks and gluons with rapidity |y| < 5 are clustered into infrared-safe jets.
After recombination we impose standard selection cuts on transverse momenta and rapidities of charged leptons and b jets, missing transverse momentum and distance between b jets according to (3.8). We require two b jets and two charged leptons in the final state, with bottom quarks in jets leading to b jets, and b jets: We have identified the renormalization scale with the factorization scale µ = µ R = µ F and have considered a fixed reference scale set to half the partonic threshold energy for ttH production according to Ref. [9]: Alternatively, we use a dynamical scale following Ref. [13] The scale uncertainty of the LO and NLO cross section is determined by variation of the renormalization and factorization scales µ R and µ F around the central value  = q,q. The third and fourth column give the integrated cross sections in fb for LO and NLO, resp. The upper and lower variations correspond to the envelope of seven scale pairs (µ R /µ 0 , µ F /µ 0 ) = (0.5, 0.5), (0.5, 1), (1, 0.5), (1, 1), (1, 2), (2, 1), (2,2). The last column provides the K factor with K = σ NLO /σ LO . and µ 0 = µ dyn for the fixed and dynamical scale choice, respectively.. While varying the renormalization scale in PDFs and matrix elements, the top-quark width remains fixed as computed at the top-quark mass according to (3.3). For the investigation of the scale dependence of the integrated cross section in Figure 3 we vary the scale µ up and down by a factor of eight for the LO and NLO integrated cross section for the three cases: 1) µ R = µ F = µ, 2) µ R = µ, µ F = µ 0 , 3) µ F = µ, µ R = µ 0 . While we show all three cases for the dynamical scale, we show only the first case for the fixed scale choice in Figure 3. For all other results, i.e. those in Table 1 and Figures 5-7, the scale uncertainties are determined from factor-two variations as follows. We compute integrated and differential cross sections at seven scale pairs, (µ R /µ 0 , µ F /µ 0 ) = (0.5, 0.5), (0.5, 1), (1, 0.5), (1, 1), (1, 2), (2, 1), (2,2). The central value corresponds to (µ R /µ 0 , µ F /µ 0 ) = (1, 1) and the error band is constructed from the envelope of these seven calculations.

Integrated cross section and scale dependence
In Table 1 we present the integrated cross sections with fixed (3.10) and dynamical scale (3.11) at the LHC at √ s = 13 TeV corresponding to the input parameters (3.1), (3.2) and (3.6) and the cuts as defined in (3.9). The results include only contributions of O α s α 5/2 for LO amplitudes and the corresponding O(α s ) QCD corrections. We neglect possible contributions to qq processes of O α 7/2 for LO amplitudes, which we determined to be about 2 per mille of the integrated cross section at LO for the setup described above. We also do not include partonic channels with incoming bottom quarks. At LO and using the   Since integrated cross sections and NLO effects are very similar, the following considerations hold true for both, the dynamical and fixed scale: The major contributions to the cross section originate from the gluon-fusion process, with about 70 % at LO while increasing at NLO to 76 %. The contribution of the quark-antiquark annihilation drops from about 30 % at LO to 19 %. At NLO the gluon-(anti)quark induced real-radiation subprocesses contribute about 5 % to the integrated cross section. The inclusion of NLO QCD corrections reduces the scale dependence from 31 % to 5 %.
We display the dependence of the integrated LO (blue) and NLO (red) cross sections on the values of the fixed and dynamical scale in

Limit of on-shell top quarks
To determine the effects of non-resonant and off-shell top-quark contributions on the integrated cross section we perform a numerical extrapolation to the zero-top-width limit, Γ t → 0. To this end we plot is the top-quark width at LO and NLO, resp., and extrapolate linearly to Γ t → 0, using a linear regression based on the computed LO and NLO integrated cross sections, as shown in Figure 4. The factor (Γ t /Γ LO/NLO t ) 2 restores the physical top-decay branching fraction. Finite-top-width effects can be extracted by comparing the results forσ LO/NLO (Γ t → 0) to σ LO/NLO (Γ LO/NLO t ). At the LHC at √ s = 13 TeV for fixed scale µ 0 = µ fix finite-top-width effects shift the LO and NLO cross section by −0.07 ± 0.01 % and −0.14 ± 0.22 %, respectively, which are within the expected order of Γ t /m t . The strong suppression of finite-top-width effects is related to the requirement of a final state with two hard b jets. Finite-width effects are much more sizable in calculations where phase-space regions allowing for associated single-top plus W-boson production are included [65]. Such calculations require massive bottom quarks to regularize collinear singularities.
Most of the displayed differential distributions were obtained using the dynamical scale µ dyn , except for Figure 5 which illustrates the effect of the fixed-scale choice on transversemomentum distributions. Where appropriate we compare NLO effects in differential distributions to those of the related process pp → e + ν e µ −ν µ bb presented in Ref. [47]. In general we find similar NLO effects for most of the distributions, but being often more distinct for pp → e + ν e µ −ν µ bb in Ref. [47].
Figures 5a and 5b display the transverse-momentum distributions of the positron and the harder b jet, resp., for the fixed scale µ 0 = µ fix . The K factor of the transversemomentum of the positron drops by about 50 % within the plotted range. In the high-p T tail the NLO predictions move outside the LO band with a scale variation of almost the same size as the LO one. The K factor of the transverse-momentum of the harder b jet exhibits the same tendency, but not as drastic as for the positron. It moves outside the LO band for p T < 60 GeV with larger scale variation as at the average p T of around 90 GeV.
In Figure 6 we collect several transverse momentum distributions obtained using the dynamical scale µ 0 = µ dyn . Figures 6a and 6c show the transverse-momentum distributions of the positron and the harder b jet, resp., to compare with the fixed-scale distributions in Figures 5a and 5b described above: They show clearly that the dynamical-scale choice improves the perturbative stability. The K factor changes only slightly (within 20 %) over the displayed range, and the NLO band lies within the LO band. The residual scale variation is at the level of 10 % at NLO.
In Figure 6b the distribution of missing transverse momentum, defined as p T,miss = p T,νe + p T,νµ , is shown. The K factor rises for p T,miss 100 GeV up to about 1.5. Figures 6c and 6d display the distribution of the transverse momentum of the harder and softer b jet, resp. While the K factor for the harder b jet exhibits a minimum at the maximum of the distribution and slightly rises towards its tail, the K factor of the softer b jet decreases by about 30 % in the plotted range.
The distribution of the transverse momentum of the b-jet pair in Fig. 14     is due to the fact that in narrow-top-width approximation the b quarks are boosted via their parent top and antitop quark, which have opposite transverse momenta resulting in a suppression of a bb system with high p T at LO. The lesser stringent kinematical constraints at NLO result in an enhancement of the cross section and thus a large K factor for high p T,bb . For the ttH production at hand a Higgs boson acquiring transverse momentum softens the kinematical constraint already at LO leading to smaller NLO corrections for high p T,bb , which can be seen in Figure 6e. Furthermore, the K factor of the distribution resembles the K factor of the missing transverse momentum due to a kinematically similar configuration, but with a stronger increase to a value of 1.8 at p T ≃ 400 GeV. Figure 6f displays the transverse-momentum distribution of the Higgs boson. The average p T of the Higgs boson is around 70 GeV. The cross section decreases more moderately with p T in the plotted range than for other transverse momentum distributions presented in Figure 6.
In Figure 7 we present further differential distributions for other types of observables: the invariant mass of the ttH and b 1 b 2 system in Figure 7a and Figure 7b, resp., the cosine of the angle between the two charged leptons in Figure 7c, the azimuthal angle in the transverse plane between them in Figure 7d, and the rapidity of the top quark and the Higgs boson in Figure 7e and Figure 7f, resp. Large NLO corrections appear in the M ttH distribution below the ttH threshold. These arise dominantly from real gluon bremsstrahlung contributions, where the emitted gluon is not recombined with the bottom quarks and thus does not contribute to M ttH . The distribution of the azimuthal angle in the transverse plane between the two charged leptons exhibits sizeable NLO effects for small angles similarly as the distribution in the cosine of the angle between the two charged leptons. The K factor varies by 40 % for these distributions. The rapidity distribution of the Higgs boson features about the typical NLO effects in the central detector region, which disappear going into forward or backward direction. The K factor of the top-quark rapidity distribution on the other hand is almost flat at the level of the NLO effects of the integrated cross section.

Checks
We have performed several comparisons and consistency checks of our calculation. We have reproduced the LO hadronic cross sections with MadGraph5_aMC@NLO [60] using the fixed scale choice. We also compared the virtual NLO contribution to the squared amplitude, 2 Re M * 0 M 1 , for the gg andūu subprocesses computed by Recola [46] for many phase-space points with MadGraph5_aMC@NLO. In Figure 8a we plot the cumulative fraction of events with a relative difference larger than ∆ between 2 Re M * 0 M 1 obtained by Recola and MadGraph5_aMC@NLO. The median of agreement is about 10 −7 for the gg subprocess and roughly 5 × 10 −9 for theūu subprocess. The agreement is worse than 10 −3 for less than 0.3 % of gg-and less than 0.04 % ofūu-subprocess events. For the gg process we checked the accuracy of virtual matrix elements in satisfying the Ward identity, by replacing the polarization vector of one of the initial-state gluons in the one-loop amplitude M 1 by its momentum normalized to its energy: ǫ µ g → p µ g /p 0 g . In Figure 8b we plot the cumulative fraction of events with 2 Re M * 0 (ǫ g )M 1 (ǫ g → p g /p 0 g )/2 Re M * 0 (ǫ g )M 1 (ǫ g ) > ∆ for virtual events obtained by Recola with a median of about 10 −9 .
With the Monte Carlo code we developed for the calculation of the process pp → e + ν e µ −ν µ bbH, which employs Recola for the evaluation of matrix elements, we reproduced the results of Ref. [47] for the closely related process pp → e + ν e µ −ν µ bb for the LHC at √ s = 8 TeV using a dynamical scale. Since the Catani-Seymour dipoles [42,44] are the same for the process without the Higgs boson, this comparison allows us to test the implementation of the subtraction formalism. As we use the same renormalization procedure as in Ref. [47], this is also verified via this comparison. In Ref. [47] the one-loop matrix elements and the Ioperator are computed in double-pole approximation for the two W-boson resonances using physical (i.e. real) W-and Z-boson masses and Γ W = Γ Z = 0. Moreover, the M H → ∞ limit is adopted, i.e. closed fermion loops involving top quarks coupled to Higgs bosons are neglected. In contrast, the new code computes the full one-loop matrix elements and the I-operator without any approximation. The largest difference between both calculations results from the use of the double-pole approximation for the virtual corrections and is of order α s Γ W /(πM W ) ∼ 0.1 %. We could reproduce the results of Ref. [47] for the integrated cross section at NLO within (2−3)σ, which is at the level of 0.5 %, thus confirming both calculations basically within the accuracy of the numerical integration. We compared our predictions for the NLO total hadronic cross section with computations of ttH production without decays of the top quarks at NLO for fixed [9] and dynamical scale [13]. To this end we performed NLO computations using the references' parameters and PDF mappings and a minimal set of cuts on the decay products of the top and antitop in our process. Since our calculation includes background processes to ttH production, a minimal set of cuts has to be maintained to ensure infrared safety: to avoid large contributions from possible collinear events with bottom quarks in forward direction (as induced by diagrams shown in Figure 1g or Figure 1h) we kept a small b-jet transverse-momentum cut (p T,b > 2 GeV) and a small b-jet distance cut (∆R bb > 0.01) to avoid singular events from gluon splitting (g → bb) as in Figure 1i. We multiply appropriate branching ratios to the results of Refs. [9,13] in the narrow-top-width approximation (NtWA), i.e. dσ NtWA = σ ttH BR t→i BRt →j , and apply the NLO matching prescription of Section 2.1.2 in Ref. [47] to our results. Thus, we find agreement of our NLO predictions with Ref. [9] and Ref. [13] within 1 %, which is of the expected order of Γ t /m t , since our calculation includes off-shell and non-resonant top-quark effects.

Conclusions
In this article we have presented the calculation of the next-to-leading-order QCD corrections to off-shell top-antitop production in association with a Higgs boson with leptonic decay of the top quarks at the LHC, including all resonant, non-resonant, and off-shell effects of top quarks as well as all interferences. For the computation of leading-and nextto-leading-order matrix elements we have utilised the recursive matrix-element generator Recola linked to the one-loop integral library Collier. The phase-space integration has been performed with an in-house multi-channel Monte-Carlo code that implements the dipole subtraction formalism. We provided integrated cross sections and several differential distributions for a 13 TeV LHC using a fixed and a dynamical scale that have both been used in the literature for the computation of NLO QCD corrections of ttH production with a stable Higgs boson and stable top quarks. We find almost the same integrated cross sections and scale dependence for both scale choices at leading order as well as next-to-leading order QCD, with a similar K factor of 1.172 and 1.176 for the dynamical and fixed scale choice, resp. However, the use of the dynamical scale instead of the fixed scale improves the perturbative stability in high-energy tails of most distributions, especially those of transverse momenta. Using the dynamical scale, we find K factors in the range 1.0−1.4 and residual scale uncertainties are the level or 10 % for distributions. For the integrated cross section an extrapolation to the zero-top-width limit has been performed, indicating non-resonant and off-shell top-quark effects below one per cent. While this effect is small, our calculation is also the first one to include NLO correction effects in the top-antitop-Higgs-boson production and the top decay processes.
Besides its phenomenological relevance, this calculation demonstrates the power of the tools Recola and Collier in performing complicated NLO calculations.